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P.I.  :  Hamid  Krim 


Our  research  effort  over  the  last  three  years,  has  primarily  focused  on  ex¬ 
ploring  of  combining  topological  information  with  geometrical  features  to  char¬ 
acterize  and  represent  3-dimensional  objects  for  classification/recognition  and 
efficient  storage  in  data  bases.  This  is  specifically  carried  out  (3-D  object  or 
graphs  of  images)  by  well  designed  mappings  from  manifolds  onto  R  and  referred 
to  as  height  functions.  While  the  topology  may  indeed  be  parsimoniously  and 
well  reflected  by  these  mappings,  much  of  the  local  geometry  is  missing.  These 
height  function  mappings  are  indeed  commonly  defined  as  equivalence  relations 
which  result  in  a  quotient  space  of  level  sets  of  the  surface,  and  which  are  eventu¬ 
ally  graphically  represented  by  the  so-called  Reeb  Graphs.  These,  as  described 
in  more  detail  later,  are  indeed  but  a  set  of  nodes  and  edges  which  bifurcate 
whenever  a  change  of  topology  takes  place.  A  natural  follow  up  investigation  is 
in  attempting,  perhaps  at  a  slight  additional  cost,  to  further  refine  such  a  rep¬ 
resentation  by  trying  to  add  on  geometrical  (more  local)  information  as  often 
required  by  problems  in  classification  and  recognition.  Towards  that  end,  and  as 
later  described,  we  introduce  a  generalized  notion  of  a  height  function^  together 
with  a  shading  function  as  well  as  a  support  function^  respectively  reflecting 
topology,  shape  and  continuity  of  the  surface  under  analysis.  In  the  process, 
we  propose  to  modify  the  existing  Reeb  graph  to  a  pseudo- weighted  graph,  as 
a  novel  compact  representation  well  adapted  to  classification  and  recognition 
problems.  The  funding  has  helped  our  group  in  completing  a  productive  and 
successful  three  years  in  publications  as  well  as  productive  on  the  publication 
front  [5,  12,  8,  11,  17,  18,  19,  20,  9,  10,  13,  14,  15,  16,  21,  22,  4,  2,  1,  1,  6,  7].  Our 
work  has  been  recognized  at  a  national  and  international  level  and  has  resulted 
in  several  invitations  to  INRJA  (Sophia- Antipolis  France),  Liapunov  Institute 
(Moscow  State  University),  and  a  3-Dimensional  Imaging  and  Computer  Vi¬ 
sion  workshop  at  University  of  Padova  (Padova, Italy),  the  Institute  of  Pure  and 
Applied  Mathematics  at  UCLA.  In  addition,  it  has  helped  us  in  securing  funds 
from  NATO  to  initiate  collaborations  between  researchers  at  Moscow  State  Uni¬ 
versity  (Dr.  Minlos)  and  Institute  for  Information  and  Transmission  Problems 
in  Moscow,  and  INRJA  at  Sophia- Antipolis  (Prance).  In  addition,  a  signifi¬ 
cant  byproduct  of  this  work  has  greatly  contributed  to  a  concerted  effort  to 
fight  crime,  in  a  project  of  Picture  ID  compression  performed  as  a  community 
outreach  effort  for  the  North  Carolina  Criminal  Justice  Information  Network 
(NCCJIN),  The  reduction  of  a  picture  ID,  thanks  to  the  hierarchical  image 
representation  inspired  by  our  graph  representations,  to  about  a  800  byte  rep¬ 
resentation  would  indeed  allow  NCCJIN  to  safeguard  their  multimillion  dollar 
investment  in  their  low  capacity  radio  network  and  manage  to  rapidly  access 
data  bases  across  the  state  and  the  nation.  Test  are  still  under  way  to  transition 
this  work  and  potentially  serve  as  an  example  for  the  police  forces  across  the 
nation  or  even  the  world.  This  work  has  been  featured  in  several  newspapers 
and  magazines  including  the  Design  Engineer,  Prism  magazine  etc. 

We  include  a  Ph.D.  dissertation  as  a  detailed  appendix  to  much  of  the  work 
funded  by  this  grant. 
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Introduction 

With  the  increasing  use  of  scanners  to  create  images,  there  is  a  rising  need  for  robust  image  denoising 
to  remove  inevitable  noise  in  the  measurements.  Even  with  high-fidelity  scanners,  the  acquired 
images  are  invariably  noisy,  and  therefore  require  filtering.  For  instance,  images  extracted  from 
voliune  data,  that  is  obtained  by  MRI  or  CT  devices,  often  contain  amotmts  of  noise  that  must 
be  removed  before  further  processing.  Removing  noise  while  preserving  the  details  is,  however, 
no  trivial  matter  since  sharp  features  are  often  blurred  and  therefore  efficient  image  denoising 
techniques  are  needed.  In  this  thesis,  we  propose  robust  variational  models  for  image  denoising 
by  solving  partial  differential  equations.  The  core  idea  behind  our  proposed  approaches  is  to  use 
geometric  insight  in  helping  construct  regularizing  functionals  and  avoiding  a  subjective  choice  of 
a  prior  in  maximum  a  posteriori  estimation.  Using  tools  from  robust  statistics  and  information 
theory,  we  show  that  we  can  extend  this  strategy  and  develop  two  gradient  descent  flows  for  image 
denoising  with  a  demonstrated  performance. 

The  major  part  of  this  thesis  is  devoted  to  a  joint  exploitation  of  geometry  and  topology  of 
objects  for  as  parsimonious  as  possible  representation  of  objects  and  its  subsequent  application 
in  object  classification  and  recognition.  The  key  idea  consists  of  capturing  geometry  along  all 
topologically  homogeneous  parts  of  an  object  by  way  of  level  curves  superimposed  on  a  Reeb  graph 
usually  extracted  by  way  of  the  object  critical  points.  A  Reeb  graph  is  a  topological  representation 
of  the  connectivity  of  a  surface  between  critical  points  which  represent  the  nodes  of  the  graph, 
and  the  edges  of  the  graph  represent  the  connected  components  of  the  surface.  Specifically,  given 
a  function  defined  on  a  surface  or  2-manifold,  a  Reeb  graph  may  be  used  to  track  its  connected 
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components  as  the  pre-image  of  the  function.  For  the  example  of  using  a  height  function,  the 
Reeb  graph  would  contain  nodes  for  each  of  the  contours  for  each  level  set  generated  by  the 
height  function.  The  resulting  skeletal  representation,  however,  is  not  rotationally  invariant  due 
to  the  rotational  non-invariance  of  the  height  function.  We  propose  a  new  methodology  called 
geodesic  shape  distribution  that  lifts  this  limitation  and  which  we  apply  to  three-dimensional  object 
matching.  The  central  idea  is  to  encode  a  3D  shape  into  a  ID  geodesic  shape  distribution.  Object 
matching  is  then  achieved  by  calculating  an  information-theoretic  measure  of  dissimilarity  between 
resulting  probability  distributions.  That  is,  the  dissimilarity  computations  are  carried  out  in  a 
low-dimensional  space  of  geodesic  shape  distributions. 

1.1  Framework  and  motivation 

1.1.1  Image  denoising 

Image  denoising  refers  to  the  process  of  recovering  an  image  contaminated  by  noise.  The  challenge 
of  the  problem  of  interest  lies  in  faithfully  recovering  the  original  image  from  the  observed  image, 
and  furthering  the  estimation  by  making  use  of  any  prior  knowledge/assumptions  about  the  noise 
process.  The  problem  of  image  denoising  has  been  addressed  using  a  munber  of  different  techniques 
including  wavelets  [49],  order  statistics  based  filters  [17],  PDE-based  algorithms  [30,76],  and  varia¬ 
tional  approaches  [26,27,6].  In  particular,  a  large  number  of  PDE-based  methods  have  particularly 
been  proposed  to  tackle  the  problem  of  image  denoising  [4,23,25]  with  a  good  preservation  of  edges. 
Much  of  the  appeal  of  PDE)-based  methods  lies  in  the  availability  of  a  vast  arsenal  of  mathematical 
tools  which  at  the  very  least  act  as  a  key  guide  in  achieving  numerical  accuracy  as  well  as  stability. 
Partial  differential  equations  or  gradient  descent  flows  are  generally  a  result  of  variational  problems 
using  the  Euler-Lagrange  principle  [36]. 

Most  medical  imaging  methods,  for  instance,  produce  full  three-dimensional  volumes.  Tradi¬ 
tionally,  the  medical  scans  are  viewed  as  a  series  of  superposed  two-dimensional  slices  of  the  full  3D 
volume,  and  are  also  perturbed  by  noise  in  the  course  of  acquisition,  transmission  or  processing. 
Figure  1.1  illustrates  a  2D  noisy  slice  of  a  3D  MR  head. 
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3D  object  MR  Image  Noisy  Image 

Figure  1.1:  Noisy  MR  image. 


1.1.2  Object  recognition 

Three-dimensional  objects  consist  of  geometric  and  topological  data,  and  their  compact  represen¬ 
tation  is  an  important  step  towards  a  variety  of  computer  vision  applications  including  indexing, 
retrieval,  and  matching  in  a  database  of  3D  models.  The  latter  will  be  the  focus  of  Chapter  6,  and 
the  motivation  behind  considering  3D  objects  is  illustrated  in  Figure  1.2. 

3D  models  do  not  depend  on  the  configiuration  of  cameras,  light  sources,  or  surrounding  objects. 
As  a  result,  they  do  not  contain  reflections,  shadows,  occlusions,  projections,  or  partial  objects, 
which  in  turn  greatly  simplifles  flnding  matches  between  objects  of  the  same  type.  For  example,  it 
is  plausible  to  expect  that  the  3D  model  of  a  Boeing747  contains  exactly  four  engines.  In  contrast, 
any  2D  image  of  this  Boeing747  may  contain  fewer  than  four  engines  (if  some  of  the  engines  are 
occluded),  or  it  may  contain  “extra  engines”  appearing  as  the  result  of  shadows. 

In  other  respects,  representing  and  processing  3D  models  is  more  complicated  than  for  sampled 
multimedia  data.  The  main  difficulty  is  that  3D  surfaces  rarely  have  simple  parameterizations. 
Since  3D  surfaces  can  have  arbitrary  topologies,  many  useful  methods  for  analyzing  other  media 
(e.g.,  Fourier  analysis)  have  no  obvious  analogs  for  3D  surface  models.  Moreover,  the  dimension¬ 
ality  is  higher,  and  this  makes  searches  for  pose  registration,  feature  correspondences,  and  model 
parameters  more  difficult. 

In  order  to  perform  3D  matching  and  to  carry  out  the  experiments,  first  we  need  a  database  of 


1.1  Framework  and  motivation 


4 


•  2D  provides  the  grayscale/color  information  in  the 
plane:  lost  of  depth  information 

•  e.g.:  2D  images  of  F-16s  and  MiG-23s  look  very 
similar,  but  in  3D  are  different 

•  3D  is  much  more  effective  for  recognition  and  dis¬ 
play 

•  3D  applications:  industry,  medicine,  search,  video 
games  and  cinema 


Figure  1.2:  Motivation  of  3D  matching. 

3D  models,  and  a  small  subset  of  oiu:  large  database  is  depicted  in  Figure  1.3.  We  collected  several 
htmdred  models  which  consist  of  military  objects,  human  body  parts,  animals  and  other  objects. 

There  are  two  major  techniques  for  3D  object  recognition:  feature-based  and  global  methods 
as  depicted  in  Figure  1.4.  Most  three-dimensional  shape  matching  techniques  proposed  in  the 
literature  of  computer  graphics,  computer  vision  and  computer-aided  design  are  based  on  geometric 
representations  which  represent  the  features  of  an  object  in  such  a  way  that  the  shape  dissimilarity 
problem  reduces  to  the  problem  of  comparing  two  such  object  representations.  Feature-based 
methods  require  that  featmres  be  extracted  and  described  before  two  objects  can  be  compared. 

An  alternative  to  feature-based  representations  is  global  methods.  The  idea  here  is  to  represent 
an  object  by  a  globed  measure  or  shape  distribution  defined  on  the  surface  of  the  object.  The  shape 
matching  problem  is  then  performed  by  computing  a  dissimilarity  measure  between  the  shape 
distributions  of  two  arbitrary  objects. 
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Figure  1.4:  3D  object  matching  diagram. 

1.1.3  Joint  exploitation  of  geometry  and  topology 

Although  topology  is  the  study  of  the  “shape”  of  curves  and  surfaces,  topology  typically  is  not 
concerned  with  the  embedding  of  that  curve  or  surface.  For  example,  topology  is  concerned  with 
the  fact  that  if  you  remove  a  point  from  a  circle,  it  becomes  a  line  segment.  This  is  true  whether 
the  circle  is  an  ellipse  or  whether  the  circle  has  knots  in  it.  In  computer  graphics,  one  cares  about 
the  embedding  and  geometry  of  a  surface.  If  one  were  asked  to  create  a  digital  representation  of 
a  coffee  cup,  no  one  would  be  happy  if  you  returned  a  model  that  looked  like  a  doughnut.  Even 
though  you  have  returned  an  object  with  the  correct  topological  shape,  the  geometrical  shape  is 
incorrect.  For  computer  graphics,  we  typically  are  not  concerned  with  purely  topological  aspects  of 
a  surface.  Thus,  the  algorithms  we  introduce,  are  founded  on  computational  topology,  but  consider 
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Figure  1.5:  Topological  equivalence  of  coffee  cup  and  doughnut. 

geometric  aspects  of  the  surface,  for  example,  the  geodesic  distance  is  a  global  geometric  measure. 
This  interplay  of  geometry  and  topology  is  inherent  in  the  discrete  nature  of  the  surfaces  used  in  the 
field  of  computer  graphics  where  great  care  is  taken  with  the  geometry  of  a  surface,  as  the  geometry 
plays  such  an  important  role  in  determining  the  appeareuice  of  a  surface.  Although  a  coffee  cup 
is  topologically  equivgJent  to  a  doughnut,  geometrically  the  shapes  differ  as  shown  in  Figure  1.5. 
And  the  difference  in  their  appearance  matters  greatly  when  the  goal  is  to  acciurately  represent  the 
appearance  of  real  world  objects.  Thus,  a  great  deal  of  work  in  computer  graphics  has  focused  on 
geometric  aspects  of  a  surface,  including  geometry  acquisition,  geometry  simplification,  geometry 
smoothing  and  geometry  compression.  However,  there  is  a  direct  relationship  between  the  topology 
and  the  geometry  of  a  surface  that  cannot  be  ignored.  Alternatively,  many  mathematicians  and 
computational  topologists  are  concerned  with  studying  purely  topological  properties  of  a  surface. 
This  thesis  takes  a  combined  approach  and  identifies  and  localizes  topological  features  within  a 
surface  by  mixing  topological  and  geometrical  approaches. 

The  connection  between  geometry  and  topology  is  given  by  a  topological  invariant  called  genus. 
The  genus  of  a  sturface  counts  how  many  “handles”  or  “hote”  the  surface  has,  and  two  surfaces 
having  the  same  genus  are  called  topologically  equivalent  Genus  is  a  global  invariant  for  the  surface 
and  its  scalar  value  provides  one  measure  of  the  complexity  of  the  surface  (e.g.,  a  genus  zero  shape 
is  much  less  complex  than  a  surface  with  higher  genus).  For  example,  the  coffee  cup  and  doughnut 
shown  in  Figure  1.5  have  the  same  genus  equal  to  one. 
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1.2  Contributions 

The  contributions  of  this  thesis  are  as  follows: 

^  Robust  and  efficient  variational  filters  for  image  denoising:  Using  the  theoretical  fun¬ 
damentals  of  robust  statistics,  a  variational  filter  referred  to  as  a  Huber  gradient  descent  flow 
is  proposed.  It  is  a  result  of  optimizing  a  Huber  functional  subject  to  some  noise  constraints, 
and  it  takes  a  hybrid  form  of  a  total  variation  diffusion  for  large  gradient  magnitudes  and 
of  a  linear  diffusion  for  small  gradient  magnitudes.  Using  the  gained  insight,  and  as  a  fur¬ 
ther  extension,  we  propose  an  information-theoretic  gradient  descent  flow  which  is  a  result 
of  minimizing  a  functional  that  is  a  hybrid  between  a  negentropy  variational  integral  and  a 
total  variation.  Illustrating  experimental  results  demonstrate  a  much  improved  performance 
of  the  approach  in  the  presence  of  Gaussian  and  heavy  tailed  noise. 

KP  A  Topological  Variational  Model  for  Image  Singularities:  Image  singularities  are  promi¬ 
nent  landmarks  and  their  detection,  recognition,  and  classification  is  a  crucial  step  in  image 
processing  and  computer  vision.  Such  singularities  carry  important  information  for  further 
operations,  such  as  image  registration,  shape  analysis,  motion  estimation,  and  object  recog¬ 
nition.  In  this  thesis,  we  propose  a  topological  gradient  descent  flow  for  image  singularities. 
The  approach  is  expressed  in  the  higher  order  variational  framework  as  a  minimizer  of  a  vari¬ 
ational  integral  involving  the  gradient  and  the  Hessian  matrix  of  the  height  function  defined 
on  a  manifold.  We  demonstrate  through  numerical  simulations  the  power  of  the  proposed 
technique  in  preserving  image  singularities. 

^  Topological  modeling  of  illuminated  surfaces  using  Reeb  graph:  We  introduce  a  reli¬ 
able  and  efficient  feature  based  object  representation  for  topological  modeling  of  three  di¬ 
mensional  illuminated  surfaces.  The  proposed  approach  encodes  an  object  into  the  Reeb 
graph  concept  from  computational  topology.  This  skeletal  structure  is  based  on  a  general¬ 
ized  height  function  in  the  light  direction  defined  on  an  illuminated  surface.  In  this  thesis, 
the  topological  properties  of  the  proposed  representation  are  analyzed  in  the  Morse  theoretic 
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framework,  and  its  close  relationship  to  the  shading  problem  is  also  highlighted.  Some  numer¬ 
ical  simulations  with  synthetic  and  real  3D  data  are  provided  to  demonstrate  the  potential 
of  object  singularities  in  topological  modeling. 

S’  Geodesic  Object  Representation  and  Recognition:  We  propose  a  shape  signature  that 
captures  the  intrinsic  geometric  structure  of  3D  objects.  The  primary  motivation  of  the 
proposed  approach  is  to  encode  a  3D  shape  into  a  one-dimensional  geodesic  distribution 
function.  This  compact  and  computationally  simple  representation  is  based  on  a  global 
geodesic  distance  defined  on  an  object  surface,  and  takes  the  form  of  a  kernel  density  estimate. 
To  gain  further  insight  into  the  geodesic  shape  distribution  and  its  practicality  in  3D  computer 
imagery,  some  numerical  experiments  are  provided  to  demonstrate  the  potential  and  the 
much  improved  performance  of  the  proposed  methodology  in  3D  object  matching.  This 
is  carried  out  using  an  information-theoretic  measure  of  dissimilarity  between  probabilistic 
shape  distributions. 

Distance  Function-based  Object  Recognition:  We  introduce  a  topological  approach  to 
object  recognition  using  a  distance  function.  Similarly  to  the  height  function  strategy  which 
consists  of  reconstructing  surface  from  its  cross-sections,  the  key  idea  behind  using  a  distance 
function  is  that  a  surface  may  be  reconstructed  from  its  intersections  with  concentric  spheres 
centered  at  the  centroid  of  the  underlying  surface.  As  the  distance  function  traverses  the 
surface  and  the  number  of  intersecting  contours  changes  for  various  distances  (i.e.  radii  of 
concentric  spheres)  from  the  barycenter  point,  the  topology  of  the  surface  varies  as  well. 

1.3  Thesis  overview 

The  organization  of  this  thesis  is  as  follows: 

□  The  Background  Chapter  contains  a  brief  review  of  essential  concepts  and  definitions  which 
we  will  refer  to  throughout  the  thesis,  and  presents  a  short  summary  of  material  relevant  to 
variationeJ  methods,  geometric  modeling  and  computational  topology. 
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□  In  Chapter  3,  we  present  a  variational  approach  to  maximum  a  posteriori  estimation  MAP 
estimation.  The  key  idea  behind  this  approach  is  to  use  geometric  insight  in  helping  construct 
regularizing  functionals  and  avoiding  a  subjective  choice  of  a  prior  in  MAP  estimation.  Using 
tools  from  robust  statistics  and  information  theory,  we  show  that  we  can  extend  this  strategy 
and  develop  two  gradient  descent  flows  for  image  denoising  with  a  demonstrated  performance. 

□  In  Chapter  4,  we  propose  a  geometric/topological  variational  model  to  preserve  degenerate 
image  singularities.  Such  singularities  carry  important  information  for  a  variety  of  image 
processing  and  computer  vision  operations,  such  as  image  registration,  shape  analysis,  object 
recognition,  etc.  The  approach  is  expressed  in  the  higher  order  variational  framework,  and  it 
is  based  on  a  variational  integral  involving  the  gradient  and  the  Hessian  matrix  of  the  height 
function  defined  on  a  manifold. 

□  Chapter  5  is  devoted  to  a  formulation  of  object  singularities  in  a  Morse  theoretic  setting. 
Then,  we  analyze  the  Reeb  graph  representation  in  the  shading  problem  framework,  and  we 
derive  some  relevant  theoretical  properties  of  the  height  function  in  the  light  direction  of  an 
illuminated  3D  object.  In  addition,  we  prove  that  such  a  height  function  is  closely  related  to 
the  generalized  bas-relief  transformation. 

□  In  Chapter  6,  we  propose  a  new  approach  for  object  matching  based  on  a  global  geodesic 
measure.  The  key  idea  behind  our  methodology  is  to  represent  an  object  by  a  probabilistic 
shape  descriptor  called  geodesic  shape  distribution  that  measures  the  global  geodesic  distance 
between  two  arbitrary  points  on  the  surface  of  an  object.  In  contrast  to  the  Euclidean  distance 
which  is  more  suitable  for  linear  spaces,  the  geodesic  distance  has  the  advantage  of  its  ability 
to  capture  the  (nonlinear)  intrinsic  geometric  structure  of  the  data.  The  geodesic  shape 
distribution  may  be  used  to  facilitate  representation,  indexing,  retrieval,  and  object  matching 
in  a  database  of  3D  models.  More  importantly,  the  geodesic  shape  distribution  provides  a  new 
way  to  look  at  the  object  matching  problem  by  understanding  the  intrinsic  geometry  of  the 
shape.  The  matching  task  therefore  becomes  a  one-dimensional  comparison  problem  between 
probability  distributions  which  is  clearly  much  simpler  than  comparing  3D  structures.  Object 
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matching  can  then  be  carried  out  by  dissimilarity  measure  calculations  between  geodesic  shape 
distributions,  and  is  in  addition  computationally  efficient  and  inexpensive. 

□  Chapter  7  is  devoted  to  the  distance  function-based  approach  to  topological  modeling  using 
Morse  theory.  Despite  the  theoretical  nature  of  the  results  presented  in  this  Chapter,  the 
key  idea  if  to  identify  and  encode  regions  of  topological  interest  of  3D  object  in  the  Morse- 
theoretic  framework.  The  main  motivation  behind  using  the  distance  function  is  its  rotation 
invariance  which  makes  it  more  adapted  to  object  recognition  than  the  Morse  height  function. 
We  prove  that  a  surface  may  be  reconstructed  from  its  intersections  with  concentric  spheres 
centered  at  the  barycenter  of  the  underlying  surface.  The  topological  changes  in  the  surface 
occur  as  we  increase  the  value  of  the  sphere  radius.  At  singular  points,  the  level  curves  of  the 
distance  function  may  split  or  merge  which  indicate  topological  changes.  We  also  show  that 
when  a  surface  is  embedded  in  a  sphere,  the  height  function  and  the  distance  function  are 
equivalent  in  a  Morse-theoretic  setting,  that  is  both  functions  have  the  same  singularities. 

□  In  the  Conclusions  Chapter,  we  summarize  the  contributions  of  this  thesis,  and  we  propose 
several  future  research  directions  that  are  directly  or  indirectly  related  to  the  work  performed 
in  this  thesis. 
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Chapter  2 - — — 

Background 

This  thesis  addresses  the  application  of  computationeil  geometry  and  topology  algorithmis  to  images 
and  three-dimensional  surfaces.  The  following  background  material  is  presented  to  provide  context 
for  this  work.  First,  a  representation  of  images  in  the  continuous  domain  is  presented  along  with 
some  basic  differential  operators  used  throughout  this  work.  Then,  a  discussion  of  the  various 
surface  representations  of  three-dimensional  objects  is  presented. 

2.1  Continuous  representation  of  images 

In  the  variational  setting,  an  image  is  usually  defined  in  the  continuous  domain  which  enjoys  a 
large  arsenal  of  analytical  tools,  and  hence  offers  a  greater  flexibility.  An  image  is  therefore  defined 
as  a  real-valued  function  «  :  fi  — >  R,  where  fl  is  a  nonempty,  bounded,  open  set  in  the  real  plane 
r2  (usually  fl  is  a  rectangle  in  R^). 

2.1.1  Differentiability  of  images 

A  pixel  location  in  Cl  is  denoted  by  «  =  {x,y),  and  the  gradient  of  u  is  denoted  by  Vu  =  {ux,Uy)^, 
where  Ux  and  Uy  are  the  first-order  partial  derivatives  with  respect  to  x  and  y  respectively.  These 
derivatives  are  illustrated  in  Figmes  2.1(b)  and  (c).  In  image  analysis,  the  image  gradient  defines  the 
orientation  of  an  edge  at  a  given  image  point.  The  gradient  orientation  or  direction  6  =  atan(uy/ux) 
gives  the  orientation  of  the  edge  as  shown  in  Figure  2.1(d).  An  edge  in  an  image  is  a  contour  across 
which  the  image  intensity  changes  abruptly.  Image  edges  are  usually  considered  to  be  discontinuities 
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of  the  image  intensity  function.  The  gradient  magnitude  l|Vu||  =  yju^  +  u^  gives  the  strength  of 
the  edge,  and  it  defines  an  edge  image  whose  pixels  are  bright  only  near  an  edge  as  depicted  in 
Figure  2.1(e).  To  detect  edges  of  any  orientation,  first  we  compute  the  gradient  of  the  image,  then 
compute  the  gradient  magnitude  at  every  pixel,  and  then  label  as  “edge  pixels”  all  pixels  whose 
gradient  magnitude  is  above  a  pre-determined  threshold. 

The  Hessian  matrix  of  an  image  u  is  defined  as  the  matrix  of  second-order  partial  derivatives 

V^u=(““ 

\  Uxy  Uyy  J 

and  its  Laplacian  is  defined  as  the  divergence  of  the  gradient  or  the  trace  of  the  Hessian  matrix 

Au  =  V  •  (Vu)  =  div(Vu)  =  Uxx  +  =  TV(V^u). 

Another  way  to  detect  edges  is  to  use  zero-crossing  of  the  Laplacian  which  crosses  zero  in  the 
neighborhood  of  an  edge,  and  this  technique  can  be  used  without  relying  on  a  threshold.  The 
Laplacian  image  is  depicted  in  Figure  2.1(f). 

The  Hessian  matrix  of  an  image  consists  of  three  terms  Uxxi  and  Uxy  The  Laplacian  ignores 
the  third  term  «md  returns  the  average  value  of  the  second  derivative  when  t£iking  every  orientation 
into  account.  While  the  Laplacian  ignores  one  of  them  and  considers  every  possible  orientation 
at  once,  the  Hessian  takes  all  three  terms  into  accoimt  and  is  orientation-dependent.  The  largest 
eigenvalue  of  the  Hessian  determines  the  orientation  along  which  the  second  derivative  is  maximal, 
while  the  smallest  eigenvalue  of  the  Hessian  returns  the  minimum  of  the  second  derivative. 

2.2  Three-dimensional  surfaces 

Surfaces  have  been  studied  by  mathematicians  for  centuries.  Topically  mathematicians  conceive 
of  surfaces  as  continuous,  for  example  a  surface  may  be  defined  as  a  continuous  function  of  two 
variables.  Each  surface  has  a  variety  of  attributes.  For  example,  the  surfaces  typically  used  in 
computer  graphics  are  oriented  2-manifolds  with  or  without  boimdary.  A  2-manifold  is  a  surface 
where  the  local  area  around  every  point  on  the  surface  is  Euclidean,  meaning,  around  each  point 
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(d)  atan{uy/ux)  (e)  ||Vtt||  (f)  A« 

Figure  2.1:  Image  differential  operators. 


the  surface  appears  to  be  nearly  flat.  The  world  we  live  on  is  £in  excellent  example  of  a  2-manifold. 
Manifolds  are  a  preferable  siurface  representation  because  the  surface  can  be  divided  into  regions 
called  charts  which  allow  2-manifolds  embedded  in  3D  to  be  flattened  into  a  two  dimensional  domain 
(through  parametrization).  Siufaces  used  in  computer  graphics  are  t3T)ically  oriented,  this  refers  to 
the  fact  that  the  surface  has  two  sides.  For  example,  a  sphere  has  two  sides,  while  a  Mobius  strip 
has  only  one  side.  Another  attribute  of  surfaces  is  whether  the  surface  is  closed  or  with  boundary. 
This  refers  to  the  number  of  open  boundary  components  of  a  surface.  For  example,  an  egg  shell  is 
closed  but  once  it  has  been  cracked  open,  it  becomes  a  surface  with  boimdaries. 

2.2.1  Image  graph 

To  apply  and  benefit  from  the  tools  of  geometry  in  image  analysis,  it  is  convenient  to  consider  the 
graph  of  an  image  u  which  is  a  surface  (2-dimensional  manifold)  M  C  defined  as  M  =  {(a:,  y,  2:)  : 
z  =  u(a:,  y)}  where  z  =  u{x,  y)  is  the  gray  level  at  position  (x,  y)  on  the  image  domain  D.  An  image 
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Figure  2.2:  A  facial  image  and  its  graph. 

and  its  surface  representation  are  depicted  in  Figure  2,2. 

Inspired  by  this  surface  representation,  the  image  denoising  problem  may  be  viewed  as  surface 
smoothing.  This  may  be  carried  out  by  minimizing  an  energy  functional  with  regularization  terms 
that  evolve  the  noisy  surface  to  the  optimal  solution  as  depicted  in  Figure  2.3, 

In  order  to  allow  partial  differentiation  and  consequently  all  the  features  of  differential  calculus 
on  M,  we  need  to  consider  a  smooth  image  u,  that  is  has  continuous  partial  derivatives  of  all  orders, 
so  that  the  manifold  M  is  differentiable.  Thus  the  study  of  this  differential  manifold  involves 
topology,  since  differentiability  implies  continuity.  A  common  way  to  smooth  an  image  u  is  to 
embed  it  into  a  family  of  images  known  as  scale-space  image.  For  example,  a  Gaussian  scale-space 
image  which  is  the  result  of  convolving  an  image  with  the  bivariate  Gaussian  density.  A  parametric 
representation  of  M  is  the  Monge  patch  defined  by  r  :  >  M  such  that  r{x^y)  =  (2;,y, u(x,  j/)). 

Note  that  the  patch  r  covers  all  M,  that  is,  r{n)  =  M,  and  it  is  regular,  that  is,  Vx  x  Vy  ^  0  or 
equivalently,  the  Jacobian  matrix  of  r  has  rank  2.  Figure  2.4  illustrates  a  parametric  representation 
of  a  surface. 

It  is  worth  pointing  out  that  the  Monge  patch  r  :  M  of  a  smooth  image  u  :  M  is  a 

diffeomorphism  because  it  is  a  smooth  bijection,  and  its  inverse  is  the  restriction  to  M  of  the 
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(c)  Noisy  surface  (d)  Filtered  surface 


Figure  2.3:  Image  denoising  as  surface  evolution. 


smooth  projection  tt  :  f2  x  M  — »  that  is  r~^  =  ttIm- 

Let  p  €  M,  then  there  exists  (a:,y)  €  fl  such  that  p  =  r{x,y).  Hence,  the  unit  normal  N  to  M 
is  given  by 

N{p)  =  N{r{x,  y))  = 

The  unit  normal  is  the  most  elementary  differential  characteristics  of  a  surface  and  determines  the 
tangent  plane  TpM  at  a  surface  point  p  =  r(x,  y).  The  tangent  plane  can  be  defined  as  the  set  of 
vectors  that  are  orthogonal  to  the  surface  imit  normal.  Since  r  is  regular,  it  follows  that  the  unit 
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Figure  2.4:  Paxametric  representation  of  a  surface. 

normal  is  well-defined  everywhere,  and  therefore  the  image  graph  M  is  orientable.  For  notational 
convenience,  the  unit  normal  can  be  viewed  as  a  mapping  S^,  called  Gauss  map  of  r  and 

is  defined  as  g{x,y)  =  N{r{x,y)),  where  =  {p  €  :  ||p||  =  1}  is  the  unit  sphere.  Note  that 

y)  denotes  the  values  of  the  unit  Gauss  map  at  p  =  r (rc,  y)  as  illustrated  in  Figure  2.5. 


Figure  2.5:  Illustration  of  the  Gauss  map. 
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2.2.2  Triangle  mesh 

In  computer  graphics,  3D  objects  are  usually  represented  as  a  triangle  mesh  M  =  (V,  T),  where 
V  =  is  the  set  of  vertices,  and  T  =  {ri,.,.,Tn}  is  the  set  of  triangles.  Triangle 

meshes  are  used  so  frequently  to  represent  surfaces  in  the  discrete  domain  that  most  computer 
graphics  hardware  is  optimized  to  render  triangles.  Example  of  triangle  meshes  are  depicted  in 
Figure  2.8.  For  triangulation,  we  use  the  barycentric  subdivision  illustrated  in  Figure  2,7.  This 
technique  consists  in  introducing  a  new  vertex  at  the  center  of  each  triangle  and  a  new  vertex  at 
the  midpoint  of  each  edge  and  drawing  edges  from  the  centroid  of  the  triangle  to  each  of  the  new 
midpoint  vertices  and  to  the  original  vertices. 

2.2.3  Scalar  volume  and  isosurface 

Another  common  discrete  surface  representation  is  an  isosurface  in  a  scalar  volume  as  pictured  in 
Figure  2.8.  A  scalar  volume  consist  of  a  regularly  sampled  3D  grid  of  scalar  values.  Scalar  volumes 
are  typically  acquired  from  veal  world  data  from  various  sources,  such  as,  magnetic  resonance 
imaging  (MRI)  and  computed  Tomography  (CT)  imaging.  These  popular  imaging  techniques  are 
used  in  a  variety  of  medical  and  scientific  applications  to  view  and  analyze  three  dimensional 
structures. 


'chapter  3 


Robust  Variational  Image  Denoising 

In  this  chapter,  we  present  a  variational  approach  to  MAP  estimation  [6,7].  The  core  idea  behind 
this  approach  is  to  use  geometric  insight  in  helping  construct  regularizing  functionals  and  avoiding 
a  subjective  choice  of  a  prior  in  MAP  estimation  [8].  Using  tools  from  robust  statistics  and  infor¬ 
mation  theory,  we  show  that  we  can  extend  this  strategy  and  develop  two  gradient  descent  flows 
for  image  denoising  with  a  demonstrated  performance  [8,9]. 

3.1  Introduction 

In  recent  years,  variational  methods  and  partial  differential  equations  (PDE)  based  methods  [57, 
61,1,65,2,74]  have  been  introduced  to  explicitly  account  for  intrinsic  geometry  to  address  a  variety 
of  problems  including  image  segmentation,  mathematical  morphology,  motion  estimation,  image 
classification,  and  image  denoising  [56,62,30,47,66,4].  The  latter  will  be  the  focus  of  the  present 
chapter.  The  problem  of  signal/image  denoising  has  been  addressed  using  a  number  of  different 
techniques  including  wavelets  [49],  order  statistics  based  filters  [17],  PDE-based  algorithms  [30, 
76],  and  variational  approaches  [26,27,6].  In  particular,  a  large  number  of  PDE-based  methods 
have  particularly  been  proposed  to  tackle  the  problem  of  image  denoising  [4, 23, 25]  with  a  good 
preservation  of  edges.  Much  of  the  appeal  of  PDEl-based  methods  lies  in  the  availability  of  a  vast 
arsenal  of  mathematical  tools  which  at  the  very  least  act  as  a  key  guide  in  achieving  numerical 
accuracy  as  well  as  stability.  Partial  diff’erential  equations  or  gradient  descent  flows  are  generally 
a  result  of  variational  problems  using  the  Euler-Lagrange  principle  [36].  One  popular  variational 


23 


3.2  Problem  statement 


24 


technique  used  in  image  denoising  is  the  total  variation  based  approach.  It  was  developed  in  [65]  to 
overcome  the  basic  limitations  of  all  smooth  regularization  algorithms,  and  a  variety  of  numerical 
methods  have  also  recently  been  developed  for  solving  total  variation  minimization  problems  [43,24]. 

In  the  next  section,  a  general  formulation  of  signal/image  denoising  problem  is  stated.  In 
Section  3.3,  we  briefly  recall  the  MAP  estimation  technique,  and  in  Section  3.4  we  formulate  the 
problem  of  MAP  estimation  in  the  calculus  of  variations  setting.  Section  3.5  is  devoted  to  a  robust 
variational  formulation  using  concepts  borrowed  from  robust  estimation,  followed  by  a  probabilistic 
interpretation  of  the  nonlinear  anisotropic  diffusion  in  the  the  order  statistics  framework.  In  Section 
3.6,  information-theoretic  variational  flows  based  on  the  differential  entropy  are  proposed.  In 
Section  3.7,  we  provide  experimental  results  to  demonstrate  a  much  improved  performance  of  the 
proposed  gradient  descent  flows  in  image  denoising.  Finally,  some  conclusions  and  discussions  are 
included  in  Section  3.8. 

3.2  Problem  statement 

In  all  real  applications,  measmrements  are  perturbed  by  noise.  In  the  course  of  acquiring,  transmit¬ 
ting  or  processing  a  digital  image  for  example,  the  noise-induced  degradation  may  be  dependent  or 
independent  of  data.  The  noise  is  usually  described  by  its  probabilistic  model,  e.g.,  Gaussian  noise 
is  diaracterized  by  two  moments.  Application-dependent,  a  degradation  often  yields  a  resulting 
signal/image  observation  model,  and  the  most  commonly  used  is  the  additive  one, 

U0=U  +  T/,  (1) 

where  the  observed  image  uo  includes  the  original  signal  u  and  the  independent  and  identically 
distributed  (i.i.d)  noise  process  tj. 

Image  denoising  refers  to  the  process  of  recovering  an  image  contaminated  by  noise  (see  Fig¬ 
ure  3.1).  The  challenge  of  the  problem  of  interest  lies  in  faithfully  recovering  the  underlying  sig¬ 
nal/image  u  from  ICO,  *tnd  furthering  the  estimation  by  making  use  of  any  prior  knowledge/ assumptions 
about  the  noise  process  t).  This  goal  is  graphically  and  succinctly  described  in  Figure  3.1. 
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Figure  3.1:  Block  diagram  of  image  denoising  process. 


3.3  MAP  estimation:  model-based  approach 


In  a  probabilistic  setting,  the  image  denoising  problem  is  usually  solved  in  a  discrete  domain,  and 
in  this  case  an  image  is  expressed  by  a  random  matrix  n  =  (u^j)  of  gray  levels.  To  account  for 
prior  probabilistic  information  we  may  have  for  u.  a  technique  of  choice  is  that  of  a  maximum  a 
posteriori  estimation.  Denoting  by  p{u)  the  prior  distribution  for  the  unknown  image  it,  the  MAP 
estimator  is  given  by 

u  =  argmax{logp(^^o|^^')  +  logp(n)},  (2) 

where  p{uo\^l)  denotes  the  conditional  probability  of  ap  given  u. 

A  general  model  for  the  prior  distribution  p(u)  is  that  of  a  Markov  random  field  (MRF)  which 
is  characterized  by  its  Gibbs  distribution  given  by  [33] 


p(a)  ==  -  exp 


where  Z  is  a  partition  function  and  A  is  a  constant  known  as  the  temperature  in  the  terminology 
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of  physical  systems.  T  is  called  the  energy  function  and  has  the  form  J^{u)  =  XlceC  where  C 

denotes  a  set  of  cliques  (i.e.  set  of  connected  pixels)  for  the  MRP,  and  Vc  is  a  potential  function 
defined  on  a  clique.  We  may  define  the  cliques  to  be  adjacent  pairs  of  horizontal  and  vertical  pixels. 
Note  that  for  large  A,  the  prior  probability  becomes  flat,  and  for  small  A,  the  prior  probability 
exhibits  sharp  modes. 

Markov  random  fields  have  been  extensively  used  in  computer  vision  particularly  for  image 
restoration,  and  it  has  been  established  that  Gibbs  distributions  and  MRP’s  axe  equivalent  (  e.g. 
see  [33]).  In  other  words,  if  a  problem  is  defined  in  terms  of  local  potentials  then  there  is  a  simple 
way  of  formulating  the  problem  in  terms  of  MRP’s.  If  the  noise  process  77  is  i.i.d.  Gaussian,  then 
we  have 

p{uo\u)  =  Kexp  )  > 

where  X  is  a  normalizing  positive  constant,  is  the  noise  variance,  and  |  •  |  stands  for  the  Euclidean 
norm  or  for  the  absolute  value  in  the  case  of  a  scalar.  Thus,  the  MAP  estimator  in  (2)  yields 


u 


=  arg  min 

u 


+  -^\u-uo 


(3) 


Image  estimation  using  MRP  priors  has  proven  to  be  a  powerful  approach  to  restoration  and 
reconstruction  of  high-quality  images.  Its  major  drawback,  besides  its  computational  load,  is  the 
difficulty  in  systematically  selecting  a  practical  and  reliable  prior  distribution.  The  Gibbs  prior 
parameter  A  is  also  of  particular  importance  since  it  controls  the  balance  of  influence  of  the  Gibbs 
prior  and  that  of  the  likelihood.  If  A  is  too  small,  the  prior  will  tend  to  have  an  over-smoothing 
effect  on  the  solution.  Conversely,  if  it  is  too  large,  the  MAP  estimator  may  be  unstable  and  it 
reduces  to  the  maximum  likelihood  solution  as  A  goes  to  infinity.  Another  difficulty  in  using  a  MAP 
estimator  is  the  non-uniqueness  of  the  solution  when  the  energy  function  T  is  not  convex. 


3.4  A  veiriational  approach  to  MAP  estimation 

Unknown  prevailing  statistics  or  underlying  signal/image/noise  models  often  make  a  “target”  de¬ 
sired  performance  quantitatively  less  well  defined.  Specifically,  it  may  be  qualitative  in  nature  (e.g.. 
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preserve  high  gradients  in  a  geometric  setting,  or  determine  a  worst  case  noise  distribution  in  a  sta¬ 
tistical  estimation  setting  with  a  number  of  interpretations),  and  may  not  necessarily  be  tractably 
assessed  by  an  objective  and  optimal  performance  measure.  The  formulation  of  such  qualitative 
goals,  is  typically  carried  out  by  way  of  adapted  functionals  which  upon  being  optimized,  achieve 
the  stated  goal,  e.g.  a  monotonically  decreasing  functional  of  gradient  modifying  a  diffusion  [61]. 
This  approach  is  the  so-called  variational  approach.  It  is  commonly  formulated  in  a  continuous 
domain  which  enjoys  a  large  arsenal  of  analytical  tools,  and  hence  oflfers  a  greater  flexibility.  An 
image  is  therefore  defined  as  a  real- valued  function  u  :  ►  R,  and  is  a  nonempty,  boimded,  open 

set  in  R^  (usually  fl  is  a  rectangle  in  R^).  Throughout,  x  =  (xi,X2)  denotes  a  pixel  location  in 
and  11*11  denotes  the  L^-norm.  While  the  ultimate  overall  objective  in  the  aforementioned  formu¬ 
lation  may  coincide  with  that  of  a  probabilistic  formulation,  namely  the  recovery  of  an  underlying 
desired  signal  u,  it  is  herein  often  implicit  and  embedded  in  an  energy  functional  to  be  optimized. 
Generally,  the  construction  of  an  energy  functional  is  based  on  some  characteristic  quantity  spec¬ 
ified  by  the  task  at  hand  (gradient  for  segmentation,  Laplacian  for  smoothing,  etc.).  This  energy 
functional  is  oftentimes  coupled  to  a  regularizing  force/energy  in  order  to  rule  out  a  great  number 
of  solutions  and  to  also  avoid  any  degenerate  solution. 

When  considering  the  signal  model  (1),  our  goal  may  be  succinctly  stated  as  one  of  estimating 
the  underlying  image  u  based  on  an  observation  uq  and/or  any  potential  knowledge  of  the  noise 
statistics  to  further  regularize  the  solution.  This  yields  the  following  fidelity-constrained  optimiza¬ 


tion  problem 


min  J^(u) 


S.t.  Ilu  —  Uop  = 


where  .F  is  a  given  functional  which  often  defines,  as  noted  above,  the  particular  emphasis  on  the 
features  of  the  achievable  solution.  In  other  words,  we  want  to  find  an  optimal  solution  that  yields 
the  smallest  value  of  the  objective  functional  among  all  solutions  that  satisfy  the  constraints.  Using 
Lagrange’s  theorem,  the  minimizer  of  (4)  is  given  by 


u  =  argmm  |.F(u)  +  -\\u-  uolfj  ,  (5) 

where  A  is  a  nonnegative  parameter  chosen  so  that  the  constraint  ||uo  —  u|p  =  is  satisfied.  In 
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preictice,  the  parameter  A  is  often  estimated  or  chosen  a  priori. 

Equations  (3)  and  (5)  show  a  close  connection  between  image  recovery  via  MAP  estimation 
and  image  recovery  via  optimization  of  variational  integrals.  One  may  in  fact  reexpress  (3)  in  an 
integral  form  similar  to  that  of  (5). 

A  critical  issue,  however,  is  the  choice  of  the  variational  integral  .F,  which  as  discussed  later,  is 
often  driven  by  geometric  arguments.  Among  the  better  known  functionals  (also  called  variational 
integrals)  in  image  denoising  are  the  Dirichlet  and  the  total  variation  integrals  defined  respectively 
as 

V{u)  =  l  f  \Vufdx  and  TV{u)  =  f  |Vu|d®, 

2  Jsi  Ju 

where  Vu  denotes  the  gradient  of  the  image  u.  The  total  variation  method  basically  consists  in 
finding  an  estimate  u  for  the  original  image  u  with  the  smallest  total  variation  among  all  the  images 
satisfying  the  noise  constraint  ||tx  — «o|P  =  where  cr  is  assumed  known.  Note  that  the  parameter 
A  controls  the  trade-off  between  noise  removal  emd  detail  preservation. 


Figure  3.2:  Total  variation. 

The  intuition  for  the  use  of  the  total  variation  integral  is  that  it  incorporates  the  fact  that 
discontinuities  are  present  in  the  original  image  u.  It  measures  the  jumps  of  u,  even  if  it  is  discon¬ 
tinuous  as  depicted  in  Figure  3.2  (courtesy  of  [75]).  The  total  variation  method  has  been  used  with 
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success  in  image  denoising,  especially  for  denoising  images  with  piecewise  constant  features  while 
preserving  the  location  of  the  edges  exactly  [23] . 

The  Dirichlet  and  total  variation  functionals  can  be  written  is  a  generalized  form  given  by 

=  [  F{\Vu\)dx,  (6) 

Jo, 

where  F  :  R"*"  — M  is  a  given  smooth  function  called  a  variational  integrand  or  Lagrangian  [36]. 
Using  (6),  we  hence  define  a  functional 

C{u)  =  + 

F{\Vu\)  +  ^\u-uo\‘^^dx,  (7) 

which  by  the  formulation  in  (5)  becomes 

u  =  arg  min  C{u) ,  (8) 

uex 

where  X  is  an  appropriate  image  space  of  smooth  functions  like  or  the  space  BV (f2)  of  image 

functions  with  bounded  variation,  or  the  Sobolev  space  Note  that  BV{Q)  is  a 

Banach  space  with  the  norm  |1«||bv  =  ll«|lLi(n)  +  TV{u),  while  H^{^)  is  a  Hilbert  space  with  the 
norm  ||«I|^i(n)  =  ||«||^  +  llVuf. 

3.4.1  Properties  of  the  optimization  problem 

A  problem  is  said  to  be  well-posed  in  the  sense  of  Hadamard  if  (i)  a  solution  of  the  problem 
exists,  (ii)  the  solution  is  imique,  (iii)  and  the  solution  is  stable,  i.e.  depends  continuously  on  the 
problem  data.  It  is  ill-posed  when  it  fails  to  satisfy  at  least  one  of  these  criteria.  To  guarantee  the 
well-posedness  of  our  minimization  problem  (8),  the  following  result  provides  some  conditions. 

Proposition  3.1  Let  the  image  space  X  be  a  reflexive  Banach  space,  and  let  T  be 

(i)  weakly  lower  semicontinuous,  i.e.  if  for  any  sequence  in  X  converging  weakly  to  u,  we  have 
F{u)  <  lim  inf k-*oo  r{u^). 

(ii)  coercive,  i.e.  T{u)  — >  oo  as  ||it||  — >  oo. 

Then  the  functional  L  is  bounded  from  below  and  possesses  a  minimizer,  i.e.  there  exists  u  €  X 
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such  that  C{u)  =  infxC.  Moreover j  ij  T  is  convex  and  A  >  0,  then  the  optimization  problem  (8) 
has  a  unique  solution,  and  therefore  it  is  stable. 

Proof:  Prom  (i)  and  (ii)  and  the  weak  lower  semicontinuity  of  the  L^-norm,  the  functional  £ 
is  weak  lower  semicontinuous,  and  coercive,  i.e.  C{u)  — >  oo  as  ||?x||  oo. 

Let  be  a  minimizing  sequence  of  £,  i.e.  C{u^)  — >  infx  £•  An  immediate  consequence  of  the 
coercivity  of  £  is  that  must  be  bounded.  As  X  is  reflexive,  thus  converges  weakly  to  w  in  X, 
i.e  u.  Thus  C{u)  <  liminfn-^oo  £(w^)  =  infjv  £.  This  proves  that  C{u)  =  infjv  £. 

It  is  easy  to  check  that  convexity  implies  weakly  lower  semicontinuity.  Thus  the  solution  of  the 
optimization  problem  (8)  exists  and  it  is  unique  because  the  L^-norm  is  strictly  convex.  The 
stability  follows  using  the  semicontinuity  of  £  and  the  fact  that  is  bounded.  ■ 


3.4.2  Numerical  solution:  gradient  descent  flows 

To  solve  the  optimization  problem  (8),  a  variety  of  iterative  methods  such  as  gradient  descent  [65], 
or  fixed  point  method  [43,24]  may  be  applied. 

The  first-order  necessary  condition  to  be  satisfied  by  any  minimizer  of  the  functional  £  given 
by  (7)  is  that  its  first  variation  5C{u\  u)  vanishes  at  u  in  direction  of  u,  that  is 


8C{u\  v)  =  3-£(t/  +  ev)  =  0, 
ae 

c=0 

and  a  solution  u  of  (9)  is  called  a  weak  extremal  of  £  [36]. 

Using  (7)  and  (9),  we  obtain  the  first  variation  6C{u]v)  (see  Appendix  A  for  a  proof) 


6C{u\  v) 


-  /.{(' 


(|V«1), 


-Vu  •  Vvj  +  A(w  —  no)v I  dx 
Vw)  +  A(u  -  uo)  I  u  dx, 


for  all  V  6  X. 


Using  the  fundamental  lemma  of  the  calculus  of  variations,  this  vanishing  first  variation  yields 
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an  Euler-Lagrange  equation  as  a  necessary  condition  to  be  satisfied  by  minimizers  of  C.  In  mathe¬ 
matical  terms,  the  Euler-Lagrange  equation  is  given  by 

-  div  +  A(u  -  uo)  =  0,  in  Q,  (11) 

where  “div”  stands  for  the  divergence  operator.  An  image  u  satisfying  (11)  is  called  an  extremal  of 

C. 

Note  that  |Vu|  is  not  differentiable  when  Vu  =  0  (e.g.  fiat  regions  in  the  image  u).  To  overcome 
the  resulting  numerical  difficulties,  we  use  the  following  slight  modification 

|Vu|,  =  v/|Vu|2  +  €, 

where  e  is  positive  sufficiently  small. 

Proposition  3.2  Let  A  =  0,  and  S  be  a  convex  set  of  an  image  space  X.  If  the  Lagrangian  F  is 
nonnegative  convex  and  of  class  ,  then  every  weak  extremal  of  C  is  a  minimizer  of  C  on  S. 

Proof:  The  convexity  of  F  yields 

F{y)  >  F{x)  +  F'{x){y  -  x),  Vx,  y  6  R+.  (12) 

By  assumption  w  is  a  weak  extremal  of  ie.  SC{u]v)  =  0  for  all  v  €  S.  This  implies  that 
F'(|Vti|)  =  0.  Therefore,  using  (12)  we  obtain 

f  F{\Vv\)dx  >  f  F{\Wu\)dx. 

Jn  Jn 

This  concludes  the  proof.  • 

By  further  constraining  A,  we  may  be  in  a  position  to  sharpen  the  properties  of  the  minimizer, 
as  given  in  the  following. 

Proposition  3.3  Let  A  =  0,  and  S  be  a  convex  set  of  an  image  space  X.  If  the  Lagrangian  F 
is  nonnegative  convex  and  of  class  such  that  F\0)  >  0,  then  the  global  minimizer  of  C  is  a 
constant  image. 
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Proof:  Using  (12),  it  follows  that  F(|Vu|)  >  F(0).  Thus  the  constant  image  is  a  minimizer  of 
jC.  Since  S  is  convex,  it  follows  that  this  minimizer  is  global.  ■ 


Proposition  3.4  Let  A  >  0,  and  S  be  a  convex  set  of  image  space  X.  If  the  Lagrangian  F  is 
nonnegative  strictly  convex  and  of  class  ,  then  an  extremal  u  of  C  is  the  unique  minimizer  of  C 
on  S. 

Proof:  Since  u  i — >  — uop  is  strictly  convex  when  A  >  0,  then  the  functional  C{u)  is  strictly 

convex  on  S,  that  is 

C{v)  >  C{u)  +  VC{u)  ■  (u  —  u). 

By  assumption  u  is  an  extremal  of  £,  thus  C{v)  >  C{u),  for  all  w  ^  ■ 

Using  the  Euler-Lagrange  variational  principle,  the  minimizer  of  (8)  may  be  interpreted  as  the 
steady  state  solution  to  the  following  nonlinear  elhptic  PDE  called  gradient  descent  flow 

ut  =  div(g(\Vu\)Vu)  —  X{u  —  uq),  in  x  R+,  (13) 

where  g{z)  =  F'{z)lz,  with  z  >  0,  and  assumed  homogeneous  Neumann  boundary  conditions.  A 
numerical  implementation  of  this  partial  differential  equation  is  discussed  in  Appendix  B. 

3.4.3  Illustrative  cases 

The  following  examples  illustrate  the  close  connection  between  optimization  problems  of  variational 
integrals  and  boundary  value  problems  for  partial  differential  equations  in  a  no  noise  constraint 
case  (i.e.  setting  A  =  0): 

(a)  Heat  equation;  Ut  =  An  is  the  gradient  descent  flow  for  the  Dirichlet  variational  integral  T>(u). 

It  is  important  to  point  out  that  the  Dirichlet  functional  tends  to  smooth  out  sharp  jumps 
because  it  controls  the  second  derivative  of  image  intensity  i.e.  its  “spatial  acceleration” ,  and 
it  diffiises  the  intensity  values  isotropically.  Figure  3.4(b)  shows  this  blurring  effect  on  a  clean 
image  depicted  in  Figure  3.4(a). 
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(b)  Perona-Malik  (PM)  equation;  It  has  been  shown  in  [77]  that  the  PM  diffusion  Uf  —  div(,g'(|Vw|)Vu) 
is  the  gradient  descent  flow  for  the  variational  integral 

/  Fri\Vu\)dx., 

Jn 

with  sample  Lagrangians  F^^{z)  =  log  [I  F  / c^)  or  F‘^{z)  —  (l  -  exp  (-z^/c^)), 
where  z  G  and  c  is  a  timing  positive'  constant.  Those  Lagrangians  are  depicted  in  Fig¬ 
ure  3.3. 


Figure  3.3:  Anisotropic  Lagrangians. 

A  minimization  of  such  functionals  encourages  the  smoothing  of  hornogenuous/small  gradient 
regions  and  the  preservation  of  edges/high  gradient  regions.  Note  that  ill-posedness  of  this 
formulation  was  addressed  in  a  number  of  papers  (e.g..  see  [77]).  A  result  of  applying  the  PM 
flow  with  F}  to  the  original  image  in  Figure  3.4(a)  is  illustrated  in  Figure  3.4(c).  It  is  worth 
noting  how  the  diffusion  takes  place  throughout  the  homogeneous  regions  and  not  across  the 
edges. 

(c)  Curvature  flow:  Ut  —  div(|^^)  corresponds  to  the  total  variation  integral. 
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While  limiting  spurious  oscillations,  TV  optimization  preserves  sharp  jumps  as  is  often  encountered 
in  “blocky”  signals/images.  Figure  3.4(d)  illustrates  the  output  of  the  TV  flow. 


(a)  Original  image  (b)  Heat  flow 


(c)  Perona-Malik  flow  (d)  TV  flow 


Figure  3.4:  Image  evolution  under  (b)  the  heat  flow,  (c)  Perona-Malik  flow,  and  (d)  total  variation 
flow. 
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3.5  Robust  variational  approach 


3.5.1  Robustness  for  unknown  statistics 

In  robust  estimation,  for  example,  a  case  where  even  the  noise  statistics  are  not  precisely  known 
[42, 49]  arises.  In  this  case,  a  reasonable  strategy  would  be  to  assume  that  the  noise  is  a  member  of 
some  set,  or  of  some  class  of  parametric  families,  and  to  pick  the  worst  case  density  {least  favorable, 
in  some  sense)  member  of  that  set,  and  obtain  the  best  signal  reconstruction  for  it.  Huber’s 
e-contaminated  normal  set  Ve  is  defined  as  [42] 

Ve  =  {{1  -  e)^  +  eH  :  HeS}, 


where  is  the  standard  normal  distribution,  S  is  the  set  of  all  probability  distributions  symmetric 
with  respect  to  the  origin  and  e  €  [0, 1]  is  the  known  fraction  of  “contamination” .  Huber  found  that 
the  least  favorable  distribution  in  Ve  which  maximizes  the  asymptotic  variance  (or,  equivalently, 
minimizes  the  Fisher  information)  is  Gaussian  in  the  center  and  Laplacian  in  the  tails.  The  transi¬ 
tion  between  the  two  depends  on  the  fraction  of  contamination  e,  i.e.,  larger  fractions  correspond 
to  smaller  switching  points  and  vice  versa. 

For  the  set  Ve  of  e-contaminated  normal  distributions,  the  least  favorable  distribution  has  a 
density  function 

fniz)  =  ((1  -  e)/y/^)exp{-pk{z)), 


where  pk  is  the  Huber  M-estimator  cost  function  (see  Figure  3.5)  given  by 

if  |z|  <  k 


Pk{z)  =  { 


k\z\  — —  otherwise. 


Here  A;  is  a  positive  constant  determined  by  the  fraction  of  contamination  e  by  the  equation 

2  (M -♦(-,))  (14) 


where  $  is  the  standard  normal  distribution  function  and  (f)  is  its  probability  density  function.  It 
is  clear  that  pk  is  a  convex  function,  quadratic  in  the  center  and  linear  in  the  tails  as  illustrated  in 
Figure  3.5. 
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Figure  3.5:  Huber  function. 


Motivated  by  the  robustness  of  the  Huber  M-filter  in  a  probabilistic  setting  and  its  resilience  to 
impulsive  noise,  we  propose  a  variational  filter  which,  when  accounting  for  these  properties,  leads 
to  the  following  energy  functional 


nk{u)=  /  pk{\Vu\)dx. 

Ju 

Note  that  the  Huber  variational  integral  is  a  hybrid  of  the  Dirichlet  variational  integral  (pfc(|  Vu|)  a 
|Vup/2  as  fc  — »  oo)  and  of  the  total  variation  integral  (p*;(|Vu|)  a  |Vii|  as  A:  — >  0).  One  may  check 
that  the  Huber  variational  integral  Tlk  ■  ►  1R+  is  well  defined,  convex,  and  coercive.  It 

follows  from  Proposition  1  that  the  minimization  problem 

u  =  arg  min  (7lfc(«)  +  ^||u  -  uo|p|  =  arg  min  /  ( Pfc(|Vu|)  +  -  uop  )  da;  (15) 

ueH^(Q)  I  2  J  ueH^(0.)jQ  \  z  / 


has  a  solution.  This  solution  is  unique  when  A  >  0. 


Proposition  3.5  The  optimization  problem  (15)  is  equivalent  to 


(02 

u  =  arg  mm  < 


\Vu\  -6\  +  -\u  -  uo| 


(16) 
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Proof:  For  z  fixed,  define  =  ^0^  +  k\z  -  d\  on  R.  It  is  is  clear  that  ^  is  convex  on 
R.  It  follows  that  ’J'  attains  its  minimum  at  9q  such  that  =  0  >  0)  that  is 

00  =  ksign(z  —  k).  Thus  we  have 


^(0o)  =  < 


kz  — 


-kz 


if  z  >  fc 
if  z  =  fe 
if  z  <  —k, 


and  therefore  pk{z)  =  argminggR  ^(0).  This  concludes  the  proof. 


Using  the  Euler-Lagrange  variational  principle,  a  Huber  gradient  descent  flow  is  obtained  as 
tt(  =  div(5fc(|Vu|)Vix)  —  A(«  —  uo),  in  fl  x  R.|.,  (17) 

where  gk  is  the  Huber  M-estimator  weight  function 

if  |z|  <  k 


9k{z. 


)  =  =  J 


,  -p-T  otherwise. 

l  ki 


For  large  A:,  this  flow  yields  an  isotropic  diffusion  (heat  equation  when  A  =  0),  and  for  small  fc, 
it  corresponds  to  total  variation  gradient  descent  flow  (curvature  flow  when  A  =  0). 

It  is  worth  pointing  out  that  in  the  case  of  no  noise  constraint  (i.e.  setting  A  =  0),  the  Huber 
gradient  descent  flow  yields  a  robust  anisotropic  diffusion  [20]  obtained  by  replacing  the  diffusion 
functions  proposed  in  [61]  with  robust  M-estimator  weight  functions  [42]. 

Recently,  we  proposed  a  smooth  Huber  variational  integral  [18]  given  by 


$(u)  =  /  (p{\S/u\)dx^ 

Jn 


where  the  Lagrangian  cp  is  defined  as 

-ct 


if  t  <  —a 
if  —a  <  t  <—b 


{t  +  a)^/3  —  ct 

-  b^)/2  +  {{a  -  bf  +  36c)/3  if -b<t<b 


—{t  -  a)^/3  +  ct 
ct 


ifb<t<a 

otherwise. 
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with  a  =  3/2,  6  =  1,  and  c  =  5/4.  Its  derivative  (p'  (also  referred  to  as  influence  function  in  robust 
statistics)  is  depicted  in  Figure  3.6.  The  Huber  influence  function,  however,  is  not  differentiable  as 
shown  in  Figure  3.6.  The  differentiability  of  the  influence  function  is  of  great  importance  since  it 
implies  the  continuity  of  its  first  derivative  which  in  turn  implies  the  continuity  of  the  confidence 
intervals  in  the  data  points.  A  more  detailed  description  of  the  smooth  Huber  gradient  descent 
flow  will  be  reported  elsewhere. 


Figure  3,6;  Huber  influence  function  and  its  smooth  version. 

3.5.2  Perona-Malik  equation:  an  estimation-theoretic  perspective 

In  a  similar  spirit  as  above,  one  may  proceed  to  justify  the  Perona-Malik  equation  from  a  specific 
statistical  model.  Assuming  an  image  u  =  (uij)  as  a  random  matrix  with  i.i.d.  elements,  the  output 
of  the  Log-Cauchy  filter  [17]  is  defined  as  a  solution  to  the  maximum  log-likelihood  estimation 
problem  for  a  Cauchy  distribution  with  dispersion  c  and  estimation  parameter  6.  In  other  words. 
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the  output  of  a  Log-Cauchy  filter  is  the  solution  to  the  following  robust  estimation  problem  [17] 

min  ^2  +  (uij  —  =  nim  ^  Fc{uij  —  ^), 

id  id 

where  the  cost  function  Fc  coincides  with  the  Lagrangian  function  which  yields  the  Perona-Malik 
equation.  Hence,  in  the  probabilistic  setting  the  Perona-Malik  flow  corresponds  to  the  Log-Cauchy 
filter.  Figure  3.7  illustrates  the  performance  of  the  Log-Cauchy  filter  in  removing  heavy-tailed 
(impulsive)  noise. 


Figure  3.7:  Log-Cauchy  filtering:  (a)  contaminated  image  with  impulsive  noise,  (b)  filtered  image. 


3.6  Information  based  functionals 

3.6.1  Information  theoretic  approach 

In  the  previous  section,  we  proposed  a  least  favorable  distribution  as  a  result  of  exercising  our 
ignorance  in  describing  that  of  an  image  gradient  within  a  domain.  Another  effective  way  is  to 
adopt  a  criterion  which  bounds  such  a  case,  namely  that  of  entropy  [32].  The  maximum  entropy 
criterion  is  indeed  an  important  principle  in  statistics  for  modeling  the  prior  probability  p(«)  of  a 
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process  u,  and  has  been  used  with  success  in  numerous  image  processing  applications.  The  term  is 
often  associated  with  qualifying  the  selection  of  a  distribution  subject  to  some  moments  constraints 
(e.g.  mean,  variance,  etc.),  that  is,  the  available  information  is  described  by  way  of  moments  of 
some  known  functions  mr(u)  with  r  =  1, . . . ,  s.  Indeed  coupling  the  finiteness  of  mr(u)  for  example 
with  the  maximum  entropy  condition  of  the  data  suggests  a  most  random  model  p{u)  with  the 
corresponding  moments  constraints  as  a  most  adapted  model  (equivalently  minimizing  negentropy 
see  e.g.  [48]). 

min  f  p{u)  log  p{u)du 
s.t.  f  p(u)du  =  1 

f  mr(u)p(u)du  =  /Xr,  7*  =  1, . . . ,  s 

Using  Lagrange’s  theorem,  the  solution  of  (18)  is  given  by 

p(u)  =  ^  exp  |-  ^  Xrmr(u) I ,  (19) 

where  Ar’s  axe  the  Lagrange  multipliers,  and  Z  is  a  partition  function.  The  resulting  model  p(u) 
given  by  (19)  may  hence  be  used  as  a  prior  in  a  MAP  estimation  formulation. 

3.6.2  Entropic  gradient  descent  flow^ 

Motivated  by  the  good  performance  of  the  maximum  entropy  principle  in  image/signal  analysis 
applications  and  inspired  by  its  rationale,  we  may  naturally  adapt  it  to  describe  the  distribution  of 
a  gradient  throughout  an  image.  Specifically,  the  large  gradients  should  coincide  with  tail  events 
of  this  distribution,  while  the  small  and  medimn  ones  representing  the  smooth  regions,  form  the 
mass  of  the  distribution.  Towards  that  end,  we  write 

n(u)=  f  ff(lVuj)dx=  f  iVullogiVuld®. 

Jn  Jci 

where  H{z)  =  z\og{z),  z>0.  Note  that  —H{z)  — >•  0  as  2:  0. 

It  follows  from  the  inequality  zlog{z)  <  z^j'l  that 

|H(u)|  <  f  IVnpd®  <  ||«ll^i(n)  <00,  Vu  €  H\n), 

J  Q 
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where  ||  •  \\H^(n)  denotes  the  F^-norm.  Thus  the  negentropy  variational  integral  H  : 
is  well  defined.  Note  also  that  the  inequality  zlog{z)  <  z^/2  implies  Ti{u)  <  T>(u)^  where  'D{u)  is 
the  Dirichlet  integral.  One  may  check  that  the  Lagrangian  H  is  strictly  convex,  and  coercive,  i.e. 
H{z)  oo  as  \z\  oo.  The  following  result  follows  from  Proposition  1. 

Proposition  3.6  Let  A  >  0.  The  minimization  problem 

u  =  arg  min  |w(w)  +  ^||u  —  uo|P  1  =  arg  min  f  (  |Vix|  log  |Vu|  + -[u  uoM  d® 

ueH^ifl)  2  J  u€H^(Ct)jQ  \  ^  / 

has  a  unique  solution  provided  that  |Vw|  >1. 

Calling  upon  the  Euler-Lagrange  variational  principle  again,  the  following  entropic  gradient 
descent  flow  results 

ut  =  div  ( ^  ^  ”  X{u  -  Wo),  in  X  E+, 

V  Jvw|  ) 

with  homogeneous  Neumann  boundary  conditions.  In  addition,  this  energy  spread  of  the  gradient 
energy  may  be  related  to  that  sought  by  the  total  variation  method,  which  in  contrast  allows  for 
additional  higher  gradients. 


Proposition  3.7  Let  u  be  an  image.  The  negentropy  variational  integral  and  the  total  variation 
satisfy  the  following  inequality 

n{u)  >  TV{u)  -  1. 

Proof:  Since  the  negentropy  H  is  a  convex  function,  the  Jensen  inequality  yields 

J  /f(|V«|)<ia;  >  |V«|da:^ 

=  h(tV{u)) 

=  Ty(ii)logTF(«), 

and  using  the  inequality  z\og{z)  >  2;  —  1  for  ^  >  0,  we  conclude  the  proof.  ■ 


Remark:  (See  Figure  3.8)  The  following  inequalities  between  Huber  variational  integral,  total 
variation  and  negentropy  integral  hold 


3.6  Information  based  functionals 


42 


(i)  If  |V«|  €  [0,e]  then  n(u)  <  TV(u) 

(ii)  If  |V«1  €  (e,cx))  then  H{u)  >  TV{u) 

(iii)  If  |Vw|  6  (e*^“^,oo)  and  fc  >  2  then  li{u)  <  TZkiu), 
where  e  is  the  Euler  number  (e  =  lim„_,oo(l  +  1/n)"  «  2.71). 


Figure  3.8:  Visual  comparison  of  some  variational  integrands 


3.6.3  Improved  entropic  gradient  descent  flow 


To  sunmaarize  and  for  a  comparison  sake,  we  show  in  Figure  3.8  the  behavior  of  the  variational 
integrands  we  have  discussed  in  this  paper.  It  can  be  readily  shown  [6]  that  a  differentiable  hybrid 
functional  between  the  negentropy  variational  integral  and  the  total  variation  may  be  defined  as 


Hrviu)  = 


H{u)  if  |V«|  <  e 

2  TV (u)  —  meas(fl)e  otherwise, 


yielding  an  improved  gradient  descent  flow.  The  quantity  meas(fl)  denotes  the  Lebesgue  measure 
of  the  image  domain  Q.  Note  that  Htv  •  K  is  well  defined,  differentiable,  convex,  and 
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Using  the  Euler-Lagrange  variational  principle,  it  follows  that  the  improved  entropic  gradient 
descent  flow  is  given  by 

Ut  =  V  •  -  A(u  -  «o),  infix]R+,  (21) 

with  homogeneous  Neumann  boundary  conditions. 


Figure  3.9:  Improved  entropic  Lagrangian. 


3.7  Experimental  results 


44 


3.7  Experimental  results 

This  section  presents  simulation  results  where  Huber,  entropic,  total  variation  and  improved  en- 
tropic  gradient  descent  flows  axe  applied  to  enhance  images  corrupted  by  Gaussian  and  Laplacian 
noise. 

The  performance  of  a  Alter  clearly  depends  on  the  Alter  type,  the  properties  of  signals/images, 
and  the  characteristics  of  the  noise.  The  choice  of  criteria  by  which  to  measure  the  performance 
of  a  Alter  presents  certain  difliculties,  and  only  gives  a  partial  picture  of  reality.  To  assess  the  per¬ 
formance  of  the  proposed  denoising  methods,  a  mean  square  error  (MSE)  between  the  Altered  and 
the  original  image  is  evaluated  and  used  as  a  quantitative  measure  of  performance  of  the  proposed 
techniques.  The  regularization  parameter  (or  Lagrange  multiplier)  A  for  the  proposed  gradient 
descent  flows  is  chosen  to  be  proportional  to  signal-to-noise  ratio  (SNR)  in  all  the  experiments. 

In  order  to  evaluate  the  performance  of  the  proposed  gradient  descent  flows  in  the  presence  of 
Gaussian  noise,  the  image  shown  in  Figure  3.10(a)  has  been  corrupted  by  Gaussian  white  noise 
with  SNR  =  4.79  db.  Figure  3.10  displays  the  results  of  Altering  the  noisy  image  shown  in  Fig¬ 
ure  3.10(b)  by  Huber  with  optimal  k  =  1.345,  entropic,  total  variation  and  improved  entropic 
gradient  descent  flows.  Qualitatively,  we  observe  that  the  proposed  techniques  are  able  to  suppress 
Gaussian  noise  while  preserving  important  features  in  the  image.  The  resulting  mean  square  error 
(MSE)  computations  are  depicted  in  Table  3.1. 

The  Laplacian  noise  is  somewhat  heavier  than  the  Gaussian  noise.  Moreover,  the  Laplace  dis¬ 
tribution  is  similar  to  Huber’s  least  favorable  distribution  [42]  at  least  in  the  tails.  To  demonstrate 
the  application  of  the  proposed  gradient  descent  flows  to  image  denoising,  qualitative  and  quan¬ 
titative  comparisons  are  performed  to  show  a  much  improved  performance  of  these  techniques. 
Figme  3.11(b)  shows  a  noisy  image  contaminated  by  Laplacian  white  noise  with  SNR  =  3.91  db. 
The  MSE’s  results  obtained  by  applying  the  proposed  techniques  to  the  noisy  image  are  shown 
in  Table  3.2.  Note  that  from  Figme  3.11  it  is  clear  that  the  improved  entropic  gradient  descent 
flow  outperforms  the  other  flows  in  removing  Laplacian  noise.  Comparison  of  these  images  clearly 
indicates  that  the  improved  entropic  gradient  descent  flow  preserves  well  the  image  structures  while 
removing  heavy  tailed  noise. 


3. 7  Experimental  results 


(a)  Original  image 


(c)  Huber 


(e)  Total  Variation 


(f)  Improved  Entropic 


Figure  3.10:  Filtering  results  for  Gaussian  noise. 
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(a)  Original  image  (b)  Noisy  image 


(e)  Total  Variation  (f)  Improved  Entropic 


Figure  3.11:  Filtering  results  for  Laplacian  noise. 
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PDE 

Mean  Square  Error  (MSE) 

SNR  =  4.79 

SNR  =  3.52  SNR  =  2.34  III 

Huber 

234.1499 

233.7337 

230.0263 

Entropic 

205.0146 

207.1040 

205.3454 

Total  Variation 

247.4875 

263.0437 

402.0660 

Improved  Entropic 

121.2550 

137.9356 

166.4490 

Table  3.1:  MSE’s  computations  for  Gaussian  noise. 


PDE 

Mean  Square  Error  (MSE) 

SNR  =  6.33 

SNR  =  3.91  SNR  =  3.05  III 

Huber 

237.7012 

244.4348 

248.4833 

Entropic 

200.5266 

211.4027 

217.3592 

Total  Variation 

138.4717 

176.1719 

213.1221 

Improved  Entropic 

104.4591 

170.2140 

208.8639 

Table  3.2:  MSE’s  computations  for  Laplacian  noise 

3.8  Discussions  and  conclusions 

In  this  chapter,  we  have  explored  a  connection  between  maximum  a  posteriori  estimation  and  the 
variational  formulation  based  on  the  minimization  of  a  given  variational  integral  subject  to  some 
noise  constraints.  A  robust  filter  called  Huber  gradient  descent  flow  was  proposed.  It  minimizes 
the  Huber  variational  integral  subject  to  some  noise  constraints.  This  filter  behaves  as  the  total 
variation  anisotropic  diffusion  for  small  gradient  magnitudes  and  as  the  isotropic  diffusion  for 
large  gradient  magnitudes.  Another  filter  called  entropic  gradient  descent  flow  derived  from  the 
maximum  entropy  principle  is  proposed.  It  minimizes  the  negentropy  variational  integral  subject 
to  some  noise  constraints.  The  proposed  gradient  descent  flows  has  been  applied  to  enhance  images 
corrupted  with  Gaussian  as  well  as  Laplacian  noise,  and  it  has  been  shown  that  these  proposed 
techniques  preserve  well  details  while  removing  noise. 


Chapter  4 


Topological  Variational  Model  for  Image 
Singularities 


In  this  Chapter,  we  propose  a  geometric/topological  variational  model  to  preserve  degenerate  image 
singularities  [16].  Such  singularities  carry  important  information  for  a  variety  of  image  processing 
and  computer  vision  operations,  such  as  image  registration,  shape  analysis,  and  object  recognition. 
The  approach  is  expressed  in  a  higher  order  variational  framework,  and  it  is  a  result  of  minimizing 
a  variational  integral  defined  in  terms  of  the  gradient  and  the  Hessian  matrix  of  the  height  function 
defined  on  a  manifold. 

4.1  Introduction 

Over  the  last  decade,  there  has  been  a  flurry  of  activity  in  applying  nonlinear  partial  differential 
equations  (PDEs)  to  image  processing  and  computer  vision.  These  approaches  have  been  proposed 
to  address  the  limitations  of  linear  scale-space  and  to  tackle  a  variety  of  imaging  applications. 
Many  of  the  these  PDE-based  techniques  are  solutions  to  variational  problems  [56,7,  8],  and  are 
determined  by  geometric  quantities.  Such  evolutions  equations  have  blossomed  in  recent  years, 
with  striking  applications  to  image  denoising,  image  segmentation,  curve  evolution,  and  motion 
estimation.  One  of  the  most  important  property  to  be  investigated  through  these  geometric  flows 
is  the  local  behavior  of  image  singularities,  and  hence  the  topology  of  the  level  sets  of  the  image. 
Depending  on  whether  the  Hessian  matrix  of  the  image  is  singular  or  not,  the  critical  points  may 
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be  divided  into  two  classes:  degenerate  and  nondegenerate  respectively  [31].  The  so-called  Morse 
theory  [3]  studies  the  properties  of  a  Morse  function  (i.e.  a  function  that  has  only  nondegenerate 
singular  points)  and  it  describes  the  topology  changes  of  the  level  sets  of  this  function  at  those 
singularities.  Regular  or  noncritical  points  do  not  affect  the  number  or  genus  of  the  components 
of  the  level  sets.  It  can  be  shown  that  Morse  functions  are  dense  and  stable  in  the  set  of  all 
smooth  functions,  that  is  the  structure  of  nondegenerate  singularities  does  not  change  under  small 
perturbations.  The  basic  ingredients  of  Morse  theory  are  Morse  lemma  and  deformation  lemma. 
The  former  states  that  in  a  neighborhood  of  a  nondegenerate  singularity,  a  function  is  reduced  to 
a  quadratic  form  in  an  appropriate  system  of  coordinates,  while  the  latter  lemma  essentially  states 
that  two  level  sets  of  a  Morse  function  are  topologically  equivalent  and  can  be  deformed  onto  one 
another  if  there  is  no  singularity  between  them. 

On  the  other  hand,  degenerate  singularities  are  unstable,  and  their  local  behavior  can  be  studied 
using  catastrophe  theory  [22]  which  deals  with  such  critical  points.  The  essential  characteristics  of 
non-Morse  function  can  be  studied  by  embedding  it  into  a  smooth  family  of  functions  controlled 
by  a  certain  parameter.  This  parameterized  function  can  be  regarded  as  a  perturbed  function, 
and  topological  changes  may  happen  when  changing  the  control  parameter  in  sudi  a  way  that  a 
degenerate  singularity  may  become  nondegenerate.  The  idea  behind  this  perturbation  method  is 
to  show  how  degenerate  singularities  behave  as  the  control  parameter  changes  and  therefore  to 
reduce  the  general  problem  to  the  nondegenerate  situation  so  that  Morse  theory  can  be  applied.  If 
we  think  of  the  control  parameter  as  time,  then  geometric  evolution  equations  can  be  regarded  as 
such  a  smooth  family  of  fimctions.  These  equations  are  flows  resulting  from  geometrical  variational 
problems,  and  are  determined  by  geometric  quantities. 

Geometry  deals  with  shape,  size  and  location  of  geometric  elements,  while  topology  is  the  con¬ 
nectivity  of  the  geometric  elements.  In  other  words,  the  geometry  of  ein  object  is  its  representation 
in  space,  while  the  topology  is  the  interconnection  of  some  geometrical  objects.  Furthermore, 
topology  is  a  global  property  of  a  space  and  does  not  depend  on  local  geometrical  structure. 

This  Chapter  is  outlined  as  follows:  the  next  section  is  concerned  with  the  problem  formulation. 
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followed  by  a  geometric  approach  to  image  singularities.  Section  4.3  describes  a  topological  charac¬ 
terization  of  singularities  in  the  Morse  theory  framework.  Furthermore,  using  the  concepts  of  height 
function  and  Gauss  map,  we  show  that  almost  all  images  are  Morse  functions.  We  also  prove  that 
images  defined  on  the  same  domain  are  topologically  equivalent.  In  Section  4.4,  a  topological  vari¬ 
ational  approach  for  image  singularities  is  proposed.  Section  4.5  is  devoted  to  experimental  results 
showing  the  methodology  proposed  in  this  Chapter.  Finally  section  4.6  presents  some  conclusions. 

4.2  Problem  formulation 

Let  u  :  Ct  G  R  be  an  image,  where  fl  is  an  open  bounded  subset  of  R^  with  Lipshitz  boundary 

(usually  Q  is  a  rectangle). 

4.2.1  Geometric  analysis  of  images 

To  apply  and  benefit  from  the  tools  of  geometry  in  image  analysis,  it  is  convenient  to  consider  the 
graph  of  an  image  u  which  is  a  surface  (2-dimension^ll  manifold)  Vi  Q  R^  defined  as  ZY  =  {{x,y,  z)  : 
z  —  u{x,y)}  where  z  =  u(x,y)  is  the  gray  level  at  position  {x,y)  on  the  image  domain  fl. 

In  order  to  allow  partial  differentiation  and  consequently  all  the  features  of  differential  calculus 
on  W,  we  need  to  consider  a  smooth  image  u,  that  is  has  continuous  partial  derivatives  of  all 
orders,  so  that  the  manifold  U  is  differentiable.  Thus  the  study  of  this  differential  manifold  involves 
topology,  since  differentiability  implies  continuity.  A  common  way  to  smooth  an  image  u  is  to 
embed  it  into  a  family  of  images  known  as  scale-space  image.  For  example,  a  Gaussian  scale-space 
image  which  is  the  result  of  convolving  an  image  with  the  bivariate  Gaussian  density.  A  parametric 
representation  of  14  is  the  Monge  patch  defined  hy  r  :  —*  14  such  that  r{x,y)  =  (®)y>^(®)y))- 

Note  that  the  patch  r  covers  all  14,  that  is,  ^(n)  =  li,  and  it  is  regular,  that  is,  r®  x  Tj,  ^  0  or 
equivalently,  the  Jacobian  matrix  of  r  has  rank  2. 

It  is  worth  pointing  out  that  the  Monge  patch  r  :  fl  — *  14  of  a  smooth  image  u  :  fl  >  R  is  a 
diffeomorphism  because  it  is  a  smooth  bijection,  and  its  inverse  is  the  restriction  to  14  of  the 
smooth  projection  tt  :  x  R  fl,  that  is  =  'k\u- 
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Let  p  eU,  then  there  exists  {x,y)  G  Q  such  that  p  =  r(x,y).  Hence,  the  unit  normal  M  to  U 
is  given  by 

»/•/  \  Krt  !  \\  1'x'^'^y  1)  /'1■\ 

Ar(p)=Ar(r(x,v))  =  ^^.^^^-— .  (1) 

Since  r  is  regular,  it  follows  that  the  unit  normal  is  well-defined  everywhere,  and  therefore  the 
image  graph  U  is  orientable.  The  unit  normal  can  be  viewed  as  a  mapping  gf  :  11  5^(1),  called 

Gauss  map  of  r  and  defined  as  g{x^y)  =  M{r{x^y)^  where  5^(1)  =  {p  :  ||p||  =  1}  is  the  unit 
sphere. 


4*2.2  Image  singularities 

Let  /  :  M  — R  be  a  real- valued  function  defined  on  a  smooth  manifold  M  C  R^.  The  function  / 
is  smooth  if  the  function  /  o  r  :  f2  R  is  smooth  (in  the  ordinary  Euclidean  sense),  where  r  is 
a  smooth  regular  patch  that  covers  all  M,  that  is,  r{Q)  =  M.  A  point  Pq  on  M  is  a  singularity 
or  critical  point  of  f  if  Po  =  r{xo,yo)i  for  some  (xo,yo)  ^  sind  the  gradient  of  /  o  r  at  (xo^yo) 
vanishes,  i.e.  V(/ o r{xo^yo))  =  0. 

A  singularity  Pq  is  nondegenerate  if  the  Hessian  matrix  V^(/  o  r(a:o5  2/o))  is  nonsingular.  Oth¬ 
erwise  this  singularity  po  is  degenerate.  The  nature  of  the  Hessian  matrix  is  directly  related  to  the 
stability  of  the  singularity.  Hence,  critical  points  can  be  divided  into  two  classes:  degenerate  and 
nondegenerate.  The  main  theory  about  nondegenerate  singularities  is  Morse  theory  that  classifies 
all  such  points.  The  degenerate  singularities  are  more  difficult  to  handle  and  catastrophe  theory 
deals  with  such  points. 


4.2.3  Image  singularities:  a  geometric  approach 

A  geometric  characterization  of  image  singularities  is  described  by  curvatures  of  the  image  graph. 
For  a  smooth  image  u,  the  Hessian  matrix  provides  information  about  the  type  of  singularities, 
and  also  characterizes  the  shape  of  the  image  graph  using  the  first  and  second  fundamental  forms 
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that  can  be  formulated  respectively  as  follows  [29] 


1  = 


and 


f  e 

m 

1  -  ( 

1 

{  m 

n  j 

'  \ 

^  ^  XX  ^  ^  xy 

'  'f*yx  ^  ‘  ‘^yy 


l  +  ul  UxUy 
V'x'^y  f  “h  Uy 


Uxx  ^xy 


\/l  +  ll^'^IP  V  Uxy  Uyy 


\/l  +  l|V«|P 


Note  that  det(I)  =  1  +  ^  0,  and  therefore  the  matrix  I  is  nonsingular. 

The  matrix  S  ~  j-ij  is  called  the  matrix  of  the  shape  operator  and  is  given  by 


5  = 


1  Sll 

S12  ' 

l_  1  / 

\  S21 

522  i 

(N 

1 

II 

gi  -  fm  gm  -  fn 


The  most  important  curvatures  in  surface  theory  are  Gaussian  and  mean  curvatures  given 
respectively  by 


K  =  det(Sr)  = 


i.n  —  m?  det(V^«) 


and 


H  =  -Tr{S)  = 

2  eg-f^ 


eg  -P  (1  +  ||Vu|P)2 
2/m  —  en  —  gi  1  „  f  Vu 


=  KiK2, 


=  2^ 


+  «2 


yi+iivuii^ 


where  Ki  and  «2  are  the  principal  curvatures  (eigenvalues  of  the  matrix  S'),  that  is 
Ki  =  H-  y/H^  -K  =  1^(5) /2  -  v/(’B:(5)/2)2  _  det(S) 


and 

K2  =  H  +  Vh^-K  =  TriS)/2  +  y/{Tr{S)/2)^  -  det(S). 

The  corresponding  eigenvectors  ei  and  62  (i.e.  Sei  =  kiCi  and  Se2  =  K2C2)  called  he  principal 
directions,  and  are  given  by  ei  =  (ki  —  S22,  S21)  and  =  (k2  “ 522)  a2i)-  These  principal  curvatures 
control  the  shape  of  the  surface  neeir  an  arbitrary  point  p  as  illustrated  in  Table  4.1.  Note  that  the 
principal  directions  may  not  be  unique.  The  signs  the  principal  curvatures  can  be  used  to  classify 
the  local  structure  of  the  image  graph. 
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Table  4.1:  Local  shape  of  a  surface. 

A  point  p  €  M  where  both  principal  curvatures  are  equal  is  called  umbilic.  Assume  that 
Ki  <  K2i  that  is,  «i  and  «2  are  the  minimal  and  maximal  principal  curvatures  with  associated 
principal  directions  ex  and  respectively. 

Definition  4.1  A  point  p  €  M  is  called  a  ridge  point  if  K2  attains  a  local  positive  maximum  in  the 
direction  of  e2,  i.e. 

«2(p)  >  0,  V«2(p)-e2  =  0  and  62  •  V^K2(p)  62  <  0. 

The  point  p  €  M  is  called  a  ravine  point  if  Ki  attains  a  local  negative  minimum  in  the  direction  of 
ex,  i.e. 

/ci(p)  <  0,  VKi(p)-ei  =  0  and  ei  •  V^Ki(p)  ei  >  0. 

4.3  Topological  analysis  of  images 

The  purpose  of  Morse  theory  is  to  explain  the  presence  and  the  stability  of  critical  points  in  terms 
of  the  topology  of  the  underlying  smooth  manifold.  Topology  is  the  property  that  determines  which 
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parts  of  the  shape  of  objects  are  connected  to  which  other  parts,  while  geometry  determines  where, 
in  a  given  coordinate  system,  each  part  is  located.  The  basic  principle  is  that  the  topology  of  a 
manifold  is  very  closely  related  to  the  critical  points  of  a  smooth  function  on  that  manifold.  Morse 
theory  studies  the  questions,  what  a  manifold  knows  about  the  critical  points  of  a  Morse  function, 
and  what  the  Morse  function  knows  about  the  manifold. 

4.3.1  Image  singularities:  a  topological  perspective 

A  typical  problem  in  mathematics  involves  attempting  to  understand  the  topology,  or  large-scale 
structme,  of  an  object  with  limited  information.  This  kind  of  problem  also  occurs  in  mathematical 
physics,  dynamic  systems  and  mechanical  engineering.  Morse  theory,  when  applied  to  smooth  man¬ 
ifolds  such  as  a  sphere  or  an  image  graph,  provides  a  general  tool  of  attacking  this  problem.  Morse 
proved  a  major  result  that  generalizes  the  straightforward  result  that  the  lowest-order  nonvanish¬ 
ing  term  in  the  Taylor  series  describes  the  local  behavior  of  a  smooth  function  of  single  variable  to 
functions  of  many  variables. 

Definition  4.2  A  smooth  function  /  :  M  R  on  a  smooth  manifold  M  is  called  a  Morse  function 
if  all  its  singularities  are  nondegenerate. 

Nondegenerate  singularities  are  isolated,  that  is,  there  cannot  be  a  sequence  of  nondegenerate 
singularities  converging  to  a  nondegenerate  singularity  p.  In  other  words,  there  is  no  other  point 
in  the  neighborhood  of  p  that  is  singular.  This  fact  follows  from  the  following  Morse’s  lemma. 

Lemma  4.3  J//  :  M  — ►  R  has  a  nondegenerate  singularity  atp^  6  M,  then  there  exists  (xq,  yo)  € 
such  that  Po  =  r{xo^yo))  o,nd  f  has  the  following  representation 

f{p)  =  /  o  yo)  ±x^±y'^  =  f{po)  ±x‘^±  y^, 

for  all  p  =  r  (x,  y)  €  M,  where  r  is  a  regular  smooth  path. 

Note  that  the  only  nondegenerate  singularities  are  the  minimum,  maximum  and  saddle  points  as 
depicted  in  Figure  4.1.  On  the  other  hand.  Figure  4.2  illustrates  a  degenerate  singularity  (cusp 


4.3  Topological  analysis  of  images 


55 


777TTTTT\ 
////// 1  1  ! 

Trm 

/  /  /  /  i 

\  \  \  \  WN 

/  /  /  /  i 

\  \  \  \ 

y'  /  /  / 

\  \  \  '^  Nk'N.'Nk 

^  ^  ^  /  t 

1— W'  .f'  ^  *  • 

' ' 

V.  v  \  \ 

N  S  \  \ 

(  f  f  /  ^ 

\  \  \  \ 

t  ;  /  /  / 

\ \ \ 

f  t  t  f  /  /  /  ^ 

WWW  \  1 

t  t  f  f  /  '// 

mm 

1 1 1  f  f // 
\t  ff/.//. 

Figure  4.1:  Nondegenerate  singular  points:  Minimum,  Maximum,  Saddle.  Bottom:  corresponding 
gradient  vector  flows. 

point).  As  we  can  see,  a  cusp  curve  is  the  intersection  of  a  surface  with  the  plane  at  the  the  cusp 
point. 

Morse  lemma  says  that  near  Pq  there  is  a  smooth  change  of  coordinates  under  which  the  resulting 
Taylor  series  of  the  Morse  function  /  near  Po  the  pure  quadratic  function. 

Theorem  4.4  Morse  functions  are  stable  and  dense  in  the  set  of  all  smooth  functions.  Equiva¬ 
lently,  any  smooth  function  can  be  converted  into  a  Morse  function  as  a  result  of  a  perturbation  as 
slight  as  desired. 

This  Morse’s  theorem  says  that  a  small,  smooth  perturbation  of  a  Morse  function  yields  another 
Morse  function.  The  density  means  that  there  is  a  Morse  function  arbitrarily  close  to  any  non-Morse 
function. 
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Figure  4.2:  Degenerate  singularity:  cusp  point/curve 
4.3.2  Almost  all  images  are  Morse  functions 

A  height  function  defined  on  a  smooth  manifold  M  is  a  real-valued  function  h  :  M  ^  R  such  that 
h(x,  y,z)  =  z  for  all  (x,  y,  z)  G  M.  Hence,  h  is  the  orthogonal  projection  with  respect  to  the  ^-axis. 
Figure  4.3  shows  a  manifold  and  the  critical  points  of  its  height  function.  These  singular  points  are 
all  nondegenerate. 

In  particular,  the  height  function  defined  on  the  graph  2/  of  an  image  u  is  given  by  h{r{x,y))  — 
u{x,  y).  Hence,  the  process  for  finding  and  classifying  the  singularities  of  the  image  u  is  the  same 
as  that  for  the  singularities  of  h  by  changing  its  local  coordinates.  Let  Po  ^  nondegenerate 
singularity  of  h  on  the  graph  U  of  the  image  u.  Then,  there  exists  (xo,yo)  €  D  such  that  Po  = 
(®o>J/OjW(a;o>i/o))*  Morse  lemma  yields 

h{x,  y,  u(x,  y))  =  h(xo,  yo,  u{xo,  yo))  ±  x^  ±  y^,  V(x,  y)  €  fi. 


or  equivalently. 


«(x,  y)  =  «(xo,  yo)  ±  x^  ±  y^,  V(x,  y)  €  D. 
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Figure  4.3:  A  3-D  object  and  the  critical  points  of  its  height  function. 

Now  we  need  to  find  out  when  a  critical  point  for  the  height  function  is  nondegenerate,  since  we 
are  interested  in  functions  whose  critical  points  are  nondegenerate.  For  this  we  use  the  concept  of 
the  Gauss  map  that  assigns  to  each  point  p  €  W  the  point  on  the  unit  sphere  that  is  parallel  to  the 
unit  normal  A/’(p).  Since  U  is  smooth,  it  follows  that  the  Gauss  map  is  smooth.  The  next  result 
links  the  height  function  to  the  Gauss  map. 

Proposition  4.5  The  height  function  on  the  graph  U  of  an  image  u  is  a  Morse  junction  if  and 
only  if  the  Gauss  map  is  a  regular  patch. 

Proof:  By  definition  the  height  function  on  the  image  graph  is  u.  Thus,  let  {x,  y)  be  a  non¬ 
degenerate  critical  point  of  the  image  u,  that  is  Vu{x,y)  =  0  and  det(V^u(x,y))  ^  0.  For  the 
sake  of  simplicity,  we  project  the  Gauss  map  g  centrally  from  the  origin  to  the  plane  z  =  1  to 
get  g{x,y)  =  {—Ux,—Uy,l)  which  is  the  unnormalized  surface  normal.  Hence,  the  cross  product 
gx^9y  =  -  det(V^tt)  7^  0,  that  is  the  Gauss  map  is  a  regular  patch.  The  converse  clearly  holds  as 
well.  ® 
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4.3.3  Topological  equivalence  of  images 

Let  u,v  :  Q.  C  R  two  smooth  images  defined  on  the  same  domain  Cl,  and  let  r  :  f2  ZY 

and  a  :  — ►  V  their  respective  Monge  patches,  where  U  and  V  are  the  image  graphs  of  u  and 
V  respectively.  Since  a  Monge  patch  of  a  smooth  image  is  a  diffeomorphism,  it  follows  that  the 
map  a  o  :  U  —y  V  is  also  a  diffeomorphism.  Hence  the  image  graphs  U  and  V  are  topologically 
equivalent,  that  is  one  can  be  transformed  into  the  other. 

4.4  Image  singularity-based  flow 

Recall  that  by  definition  a  degenerate  singular  point  satisfies  the  conditions 

Vu  =  0  and  det(V^u)  =  0. 

Hence,  we  might  formulate  the  degeneracy  of  a  singularity  in  the  calculus  of  variations  framework 
[36]  as  an  optimization  problem  given  by 

mini  f  |(1  —  €)|Vul  +  e] det(V^tt)|} d®,  €€[0,1]  (2) 

^  2  Jfi 

Taking  the  first  variation  and  using  the  identity  ■^det{M)  =  det(M)(M)”^  for  a  non-singular 
matrix  M,  we  derive  the  Euler-Lagrange  equation  which  can  be  solved  using  the  following  gradient 
descent  flow  given  by 

/n  \  ,  /  det(V2u)  .  A 

ixt  =  (1  -  €)V  •  j  -  €  A  •  ’ 

where  A  =  is  a  second  order  differential  operator.  Here  we  assume  homogeneous 

Dirichlet  and  Neumann  boundary  conditions  u  =  =  0,  where  Uj/  denotes  the  directional  derivative 

of  u  in  the  direction  of  the  unit  normal  v. 

4.5  Numerical  simulations 

We  demonstrate  the  performance  of  the  proposed  technique  by  applying  it  to  real  images.  A  test 
image  of  a  fighter  jet  F16  is  shown  in  Figure  4.4(a).  In  order  to  extract  the  degenerate  critical 
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points,  we  use  a  zero  crossing  technique  to  find  the  locations  where  both  the  gradient  and  the  Hessian 
determinant  of  the  image  vanish.  Figure  4.4(b)  shows  the  extracted  330  degenerate  singularities 
that  are  indicated  by  the  sign  “+”  in  the  image.  Figure  4.4(c)  and  4.4(d)  show  the  results  of 
applying  the  heat  flow  and  the  proposed  flow  with  their  corresponding  degenerate  singularities.  As 
we  can  see,  the  proposed  flow  preserves  more  degenerate  critical  points  (178  landmarks),  while  for 
the  heat  flow  the  number  of  landmarks  drops  rapidly  to  12.  In  Figure  4.5,  we  plot  the  number  of 
degenerate  singular  points  versus  the  iteration  number  for  the  heat  flow  and  the  proposed  technique, 
and  it  can  be  seen  clearly  that  the  proposed  approach  performs  the  best. 


(a)  Original  image  (b)  330  Landmarks 


(c)  Heat  flow  (12  Landmarks)  (d)  Proposed  flow  (178  Landmarks) 


Figure  4.4:  Evolution  of  image  singularities. 
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Figure  4.5:  Number  of  degenerate  singularities  vs.  iteration  number. 

4.6  Conclusions 

We  have  presented  an  image  singularity-based  flow  expressed  in  the  calculus  of  variations  frame¬ 
work.  This  approach  is  based  on  a  variational  functional  that  incorporates  both  differential  ge¬ 
ometry  and  topology.  The  proposed  model  provides  motivation  for  further  investigation,  and  om 
future  efforts  will  be  focused  in  applying  this  approach  to  a  variety  of  image  processing  and  com¬ 
puter  vision  tasks  such  as  image  registration,  motion  estimation,  shape  from  shading,  and  object 
recognition. 


Chapter  5 


Topological  Modeling  of  Illuminated  Surfaces 

In  this  Chapter,  we  give  a  formulation  of  object  singularities  in  a  Morse  theoretic  setting.  Then, 
we  analyze  the  Reeb  graph  representation  in  the  shading  problem  framework,  and  we  derive  some 
relevant  theoretical  properties  of  the  height  function  in  the  Ught  direction  of  an  illuminated  3D 
object  [18].  In  addition,  we  prove  that  such  a  height  function  is  closely  related  to  the  generalized 
bas-relief  transformation. 


5.1  Introduction 

The  Reeb  graph  concept  is  an  efficient  representation  that  captures  the  topological  properties  of 
three-dimensional  (3D)  objects.  The  vertices  of  the  Reeb  graph  axe  the  singular  points  of  a  function 
defined  on  the  underlying  surface  of  a  3D  object  [67,31].  The  vertices  of  the  Reeb  graph  axe  the 
singular  points  of  a  function  defined  on  the  underlying  object  [67,68,31].  These  singularities  are 
prominent  landmarks  and  their  detection,  recognition,  and  classification  is  a  crucial  step  in  image 
processing  and  computer  vision  [31].  Such  singularities  carry  important  information  for  further 
operations,  such  as  image  registration,  shape  analysis,  motion  estimation,  object  recognition,  and 
surface/image  evolution  [52, 16].  Depending  on  whether  the  Hessian  matrix  of  a  function  defined 
on  the  underlying  object  is  singular  or  not,  the  critical  points  may  be  divided  into  two  classes: 
degenerate  and  nondegenerate  respectively  [55,31].  The  so-called  Morse  theory  [55,31]  studies  the 
properties  of  a  Morse  function  (i.e.  a  function  that  has  only  nondegenerate  singular  points)  and 
it  describes  the  topology  changes  of  the  level  sets  of  this  function  at  those  singularities.  Regular 
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or  noncritical  points  do  not  affect  the  number  or  genus  of  the  components  of  the  level  sets.  It 
can  be  shown  that  Morse  functions  are  dense  and  stable  in  the  set  of  all  smooth  functions,  that  is 
the  structure  of  nondegenerate  singularities  does  not  change  under  small  perturbations.  The  basic 
ingredients  of  Morse  theory  are  Morse  lemma  and  deformation  lemma  [55,31].  The  former  states 
that  in  a  neighborhood  of  a  nondegenerate  singularity,  a  function  is  reduced  to  a  quadratic  form  in 
an  appropriate  system  of  coordinates,  while  the  latter  lemma  essentially  states  that  two  level  sets 
of  a  Morse  function  are  topologically  equivalent  and  can  be  deformed  onto  one  another  if  there  is 
no  singularity  between  them.  The  Reeb  graph  representation  is  a  result  of  extracting  and  encoding 
the  singular  points  of  a  Morse  function  defined  on  a  3D  object. 

The  rest  of  this  Chapter  is  organized  as  follows.  The  next  section  is  devoted  to  the  topological 
characterization  of  object  singularities  with  an  emphasis  on  Morse  theory  and  its  implications, 
followed  by  a  mathematical  description  of  the  Reeb  graph  representation  for  3D  data.  In  Section 
5.3,  we  establish  a  link  between  the  shading  problem  and  the  height  function  in  the  light  direction. 
Then,  we  derive  some  relevant  properties  of  this  height  function,  and  we  show  its  relation  to  the 
generalized  bas-relief  transformation.  And  finally  in  Section  5.4,  we  provide  numerical  simulations 
to  show  the  application  and  the  power  of  object  singularities  in  topological  modeling  through  the 
Reeb  graph  representation. 


5.2  Reeb  graph  representation 

5.2.1  Morse  theory  and  singularities 

Morse  theory  explains  the  presence  and  the  stability  of  critical  points  in  terms  of  the  topology 
of  the  underlying  smooth  manifold.  Topology  is  the  property  that  determines  which  parts  of  the 
shape  of  objects  are  connected  to  which  other  parts  [38],  while  geometry  determines  where,  in  a 
given  coordinate  system,  each  part  is  located  [29].  The  basic  principle  is  that  the  topology  of  a 
manifold  is  very  closely  related  to  the  critical  points  of  a  smooth  function  on  that  manifold  [55]. 

Let  ^  :  M  R  be  a  real- valued  function  defined  on  a  smooth  manifold  M  C  R^.  The  function 
<p  is  smooth  if  the  composition  function  ^  o  r  :  n  — ►  R  is  smooth  (in  the  ordinary  Euclidean  sense), 
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where  r  is  a  smooth  regular  parametrization  of  M  (i.e.  r  :  ^  E^).  A  point  Pq  on  M  is  a 

singularity  or  critical  point  of  v?  if  Po  =  for  some  (xo,yo)  ^  ^^5  and  the  gradient  of  ip  or 

at  (xo,yo)  vanishes,  i.e.  V(v?  o  r(xo,yo))  =0. 

A  singular  point  Po  is  nondegenerate  if  the  Hessian  matrix  o  r(xo,yo))  Is  nonsingular 

[55,21,31]. 

Definition  5.1  A  smooth  function  99  :  M  — >  R  on  a  smooth  manifold  M  is  called  a  Morse  function 
if  all  its  singular  points  are  nondegenerate. 

Nondegenerate  singularities  are  isolated,  that  is,  there  cannot  be  a  sequence  of  nondegenerate 
singularities  converging  to  a  nondegenerate  singularity  p.  This  fact  follows  from  a  Morse  lemma 
which  says  that  near  po  there  is  a  smooth  change  of  coordinates  under  which  the  resulting  Taylor 
series  of  the  Morse  function  h  near  Pq  is  a  pure  quadratic  function.  Note  that  the  only  nondegenerate 
singularities  are  the  minimum,  maximum  and  saddle  points. 

Another  important  result  is  Morse  theorem  which  says  that  a  small,  smooth  perturbation  of  a 
Morse  function  yields  another  Morse  function.  The  density  means  that  there  is  a  Morse  function 
arbitrarily  close  to  any  non-Morse  function. 

5.2.2  Height  function 

A  height  function  on  a  smooth  manifold  M  is  a  real-valued  function  h  such  that  h{x^  y,  z)  = 

z  for  all  (x,y,  z)  G  M.  Hence,  h  is  an  orthogonal  projection  with  respect  to  the  2:-axis.  Figure  5.1 
shows  a  2D  manifold  (a  double  torus)  and  the  critical  points  of  its  height  function.  These  singular 
points  are  all  nondegenerate. 

5.2.3  Reeb  graph 

An  interesting  concept  related  to  Morse  theory  and  very  useful  to  analyze  a  surface  topology  is  the 
Reeb  graph.  The  latter  is  defined  as  a  quotient  space  M/^  with  the  equivalence  relation  given  by 
p  q  if  and  only  if  h{p)  =  h{q)  and  p,q  belong  to  the  same  connected  component  of  h~^{h{p)). 
In  particular,  each  connected  component  is  represented  by  a  point  in  the  Reeb  graph  as  illustrated 
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Figure  5.1:  A  3-D  object  and  the  critical  points  of  its  height  function. 

in  Figure  5.2.  The  left  figure  shows  a  torus  with  the  critical  points  of  its  height  function  (Morse 
function).  The  figure  in  the  middle  illustrates  the  geometric  features  of  the  torus  represented  by 
cross-sections  along  its  height.  The  right  figure  shows  the  topological  features  represented  by  the 
Reeb  graph.  By  taking  an  appropriate  number  of  cross-sections  and  smooth  interpolation  between, 
Shinagawa  et  aL  [67,  68]  proposed  a  Reeb  graph  based  approach  or  so-called  homotopy  model 
for  object  reconstruction.  The  Reeb  graph  is  a  topological  representation  of  an  object  (skeletal 
structure),  and  has  the  advantage  to  be  stored  or  transmitted  with  a  much  smaller  amount  of  data. 

Mathematically,  a  quotient  space  M/~  =  {[p]  :  p  E  M}  is  the  set  of  equivalence  classes  of  the 
relation  where  [p]  =  {q  G  M  :  q  ^  p}  is  the  equivalence  class  of  p  G  M.  Intuitively,  M/^  is  a 
space  created  by  taking  the  space  M  and  gluing  p  to  any  q  that  satisfies  q  p-  The  classes  [p]  are 
the  connected  components  for  the  Reeb  graph,  and  being  in  the  same  component  is  an  equivalence 
relation: 

g  ~  p  h{q)  =  h{p)  and  p,q  €  ConComp{/i“^(/i(p))}, 

where  ConComp{*}  denotes  the  connected  component.  In  the  Reeb  graph  representation,  each 
connected  component  of  a  contour  (i.e.  h~^{z)  where  z  =  h{x,y,z))  corresponds  to  a  point. 
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Figure  5.2:  Reeb  graph  representation  of  a  torus.  The  dot  •  denotes  a  critical  point  of  the  height 
function. 

5.3  Shading  problem  and  height  function 

5.3.1  Shading  function 

Shadows  provide  perceptually  important  information  about  the  shape  of  an  illuminated  simface. 
Shadows  occur  when  objects  totally  or  partially  occlude  direct  light  rays  from  a  light  source.  A 
shadow  can  be  divided  into  two  classes:  self-shadow  (or  attached  shadow)  and  cast  shadow  [41]. 
The  former  is  a  portion  of  a  surface  not  illuminated  by  light  rays  (i.e.  facing  away  from  the  light 
source),  while  the  latter  is  the  area  projected  by  the  object  in  the  direction  of  light  source  on  the 
surface  (extrinsic  cast  shadow)  or  projected  on  the  smface  itself  (intrinsic  cast  shadow).  Note  that 
a  convex  object  such  as  an  egg-shell  does  not  cast  shadow  on  itself,  that  is  there  is  no  intrinsic  cast 
shadow.  Both  shadows  provide  perceptually  important  information  about  the  shape  of  the  surface 
or  object.  These  two  types  of  shadows  are  illustrated  in  Figure  5.3(a) 

Let  M  be  a  Lambertian  surface  with  unit  constant  albedo,  that  is,  it  reflects  light  equally  in  all 
directions.  Assuming  distant  point  source  illumination  in  the  direction  L  =  (^i,^2>^3)-  The  imit 
vector  L  e  (unit  sphere)  is  pointing  towards  the  light  source.  We  define  the  shading  function 
<T  :  M  R  as  the  inner  product  cr(p)  =  N(p)  •  L,  where  N  is  the  surface  unit  normal  (see 
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Figure  5.3:  (a)  Self-shadow  and  cast  shadow,  (b)  illumination  of  a  Lambertian  siurface. 


Figure  5.3(b)).  It  is  important  to  note  that  the  shading  function  has  a  similar  definition  to  that  of 
the  surface  curvatures  such  as  the  Gaussian  curvature  :  M  R. 

The  zero-level  set  C  =  cr“^(0)  of  the  shading  function  is  called  the  horizon  curve,  that  is, 
the  set  of  points  where  the  light  direction  is  orthogonal  to  the  surface.  ’The  set  of  points  in  the 
surface  not  reached  by  the  light  rays  (the  part  of  the  smface  that  is  not  illuminated),  that  is, 
E  =  {p  €  M  :  a(p)  <  0},  is  called  self-shadow.  It  is  known  that  if  the  surface  M  is  convex 
and  illuminated  in  all  directions,  then  the  self-shadow  is  a  coimected  set.  The  converse  problem, 
however,  does  not  always  hold.  Ghomi  [34, 35]  recently  studied  the  question:  does  connectedness 
of  the  self-shadows  imply  convexity  of  the  surface?.  Using  Morse  theory,  Ghomi  proved  that  the 
answer  is  yes  provided  that  each  self-shadow  is  simply  connected. 

For  almost  all  L  €  5^  the  horizon  curve  is  a  regular  curve.  It  can  also  be  shown  that  if  the 
horizon  curve  is  connected,  then  it  coincides  with  the  boundary  of  the  self-shadow  [35]. 

In  particular,  if  we  consider  the  surface  M  as  a  graph  of  an  image  or  function  u  :  Q  C.  R^  —*  R, 
that  is,  M  is  a  parametric  surface  given  by  a  Monge  patch  r  :  M  — R,  then  for  notational 
convenience,  we  may  abbreviate  a  o  r  to  <t,  so  that  the  shading  function  or  shading  image  (i.e.  the 
image  of  the  Lambertian  sm-face)  may  be  defined  as  ct  :  — >  R  such  that  cr(x,  y)  =  N{x,  y)  •  L  = 
{—Uxi\  —  Uy£2,  +  “v  for  siU  {x,y)  e  11.  In  other  words,  a{x,y)  will  denote  the  values 

of  the  shading  function  at  r{x,  y).  The  shading  image  depends  on  the  illumination,  the  properties 
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of  the  underlying  surface  and  its  orientation.  It  can  be  expressed  as  the  single  a  between  the  unit 
surface  normal  N  and  the  direction  of  the  light  source  L. 

The  self-shadow  image  E  =  {(x,y)  G  :  a{x,y)  <  0}  is  defined  as  the  orthographic  projection 
of  the  self-shadow  curve  onto  the  image  plane,  that  is,  the  {x,  y)-plane. 

The  local  properties  of  the  shading  image  are  better  described  with  respect  to  its  local  Gauge 
frame  defined  in  terms  of  the  normal  and  the  tangent  to  the  level  sets  of  the  shading  image 
as 


V  = 


Va 


and  ^  = 


(Vor)^ 


|Va|  |Va| 

Figure  5.4  shows  the  normal  and  tangential  vector  fields  of  the  shading  image  for  a  synthetic  vase 
model. 


Figure  5.4:  (a)  Shading  image,  (b)  gradient  vector  field,  (c)  shading  normal  field,  and  (d)  shading 
tangent  field. 
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5.3.2  Height  function  in  the  direction  of  light 


The  height  function  in  the  light  direction  L  of  an  illuminated  surface  M  is  given  by  h{p)  =  p  -  L, 
for  all  p  €  M. 

In  particular,  if  M  is  a  parametric  surface  given  by  r  :  >  M  such  that  r{x,  y)  —  {x,  y,  u{x,  y)) 

where  u  is  a  given  image  or  a  real-valued  function,  then  the  height  function  on  M  in  the  direction 
of  light  L  is  the  composition  of  r  with  the  orthogonal  projection  to  the  line  spanned  by  L,  and  is 
defined  as  h{x,  y)  =  r{x,  y)  •  L  =  i\x  -H  l^y  -t-  (.zu{x,  y). 

Denote  by  H  the  graph  of  the  function  h,  and  let  qeH.li  is  easy  to  verify  that  q  =  Zp,  where 
Z  is  the  matrix  given  by 


Z  = 


1 

0 


0  0 
1  0 
^2  h  ) 


Proposition  5.2  Assume  that  £3  ^  0,  and  let 
then  have 


N  = 


/l  +  l|V«f 

Vl  +  l|Vh|P 


N  be  the  unit  surface  normal  the  surface  H,  we 


\  1/2 

J  det{Z)Z-'^N. 


Proof:  By  definition,  the  unit  normal  to  the  surface  H  is  given  by 


N 


Tx  X  T" y 


(\\rx  X  T-y||\ 

vFx  xfyiiy 


£3  0  -£i  ^ 


V 


0  £3  -£2 
001 


Vx  X  r 


V 


\rX  X  Ty 


where  r  and  r  are  the  Monge  patches  for  M  and  H  respectively. 
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Proposition  5.3  Let  (p^q)  G  M  x  W.  The  point  p  is  a  on  the  self-shadow  of  M  in  the  light 
direction  L  if  and  only  if  the  point  q  is  on  the  self-shadow  ofTL  in  the  light  direction  L  =  ZL. 


Proof:  The  result  follows  using  Proposition  5.2  and  the  following  relation 


1  +  IIV^IP 
i  +  lW 


2\1/2 


det(2)(iV-L). 


Thus,  iV  •  L  <  0  if  and  only  if  iV  •  L  <  0. 


It  is  interesting  to  point  out  that  the  height  function  in  the  direction  of  an  arbitrary  vector 
(a,6,c)  with  c  >  0  is  also  referred  to  as  the  generalized  bas-relief  transformation  proposed  by 
Belhumeur  et  al.  [5]  who  have  shown  that  there  exists  an  ambiguity  in  determining  the  structure  of 
the  underlying  surface  since  both  the  surface  and  its  generalized  bas-relief  surface  produce  identical 
set  of  images  under  arbitrary  illumination,  and  therefore  they  are  indistinguishable  for  recognition 
purposes  [5].  The  above  propositions  prove  that  the  generalized  bas-relief  transformation  is  not 
the  only  one  that  produces  the  same  set  of  images  under  arbitrary  illumination  as  suggested  in  [5], 
and  that  the  height  function  in  the  light  direction  on  an  illuminated  surface  is  another  alternative 
representation  to  further  understand  and  investigate  the  shading  problem. 

5.3.3  Singularities  of  the  shading  function 

The  shading  function  of  a  manifold  M  is  defined  as  cr(p)  =  N{p)  •  L,  where  N  is  the  surface  unit 
normal,  and  L  =  is  a  unit  vector  representing  the  light  source  direction.  By  a  rigid 

motion,  we  may  move  M  tangent  to  the  (a;,  j/)-plane  at  p  =  0,  so  that  M  is  locally  a  graph  given 
hy  z  =  u{x,y)  with  u(0,0)  =  Ua;(0,0)  =  %(0,0)  =  0.  Thus,  the  unit  normal  can  be  written  as 
N  =  {-Ux,  -Uy,  l)/yjl  +  ul  +  u^.  To  study  the  singularities  of  the  shading  function  and  therefore 
the  singularities  of  of  the  Gauss  map,  it  is  more  easy  to  simplify  the  expression  of  the  unit  normal 
by  projecting  centrally  from  the  origin  to  the  plane  z  =  1  to  get  1)?  then  project 

to  the  (x,  2/)-plane  to  get  a  simplified  mapping  referred  to  as  the  modified  gauss  map  N  defined  as 
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JV(x,  y)  =  {’-Ux’i  "■%)•  Recall  that  a  central  projection  of  a  point  p  =  (x,  y,  2^)  €  M  onto  the  plane 
z  is  equal  to  {x/z^  yjz.,  1). 

The  modified  gauss  map  is  singular  when  its  Jacobian  matrix  has  rank  less  than 

2,  that  is,  when  det(V^u)  =  0.  On  the  other  hand,  assuming  that  the  third  component  of  the 
light  direction  is  positive  (i.e.  ^3  >  0),  the  shading  function  is  equivalent  a  shading  function  with 
light  direction  where  i\  and  ^2  are  arbitrary  and  ^3  =  1.  Hence,  the  shading  function 

simplifies  to  (j(x,  y)  =  iV(x,  y)  •  L,  where  L  =  (^1, 4)-  The  shading  function  is  still  denoted  a{x,  y) 
for  convenience. 

Proposition  6.4  Then  the  gradient  of  the  shading  function  is  given  by 

Vo-(x,  y)  =  -(V^«(x,  y))L. 

A  critical  point  of  the  shading  function  satisfies  (V^u)L  =  0.  At  regular  points  of  the  shading 
function,  the  horizon  curve  is  locally  a  smooth  curve. 

Proof:  Taking  the  gradient  of  the  shading  function,  we  get 

Va{x,y)  =  (-Uxxh  -  Uxy£2,  -Uxy^i  -  Uyy£2)  =  -(V^«(x,y))L. 

Hence,  Vcr(x,  y)  =  0  if  and  only  if  {V‘^u)L  =  0.  At  a  regular  point,  we  have  {V%)L  ^  0,  and  using 
the  implicit  function  theorem,  it  follows  that  the  horizon  curve  is  smooth  in  the  neighborhood  of  a 
regular  point.  ■ 


The  above  result  shows  that  the  local  orientation  of  the  shading  fimction  is  related  to  the  local 
curvature  of  the  underlying  surface.  Furthermore,  in  the  orthonormal  coordinate  frame  {61,62}, 
where  61  and  62  are  the  principal  directions,  the  Hessian  matrix  can  be  expressed  as 


f  til  0 
Y  0  K2 


1 


where  /ci  and  «2  are  the  principal  cmvatures.  Hence  the  grjidient  of  the  shading  function  be¬ 
comes  Vcr  =  — —  Ki^i.  At  a  singular  point  of  the  shading  function,  we  have  atan(K2/«i)  = 
—  atan(£i/£2),  that  is  the  shape  index  of  the  surface  depends  only  on  the  light  coordinates. 
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The  next  result  gives  a  necessary  and  sufficient  condition  for  a  nondegenerate  singular  point  of 
the  shading  function. 

Proposition  5.5  The  Hessian  matrix  of  the  shading  function  is  given  by 

o  /  \  i  ~'^xxx^l  '^xxy^2  '^'^xxy^l  '^xyy^2 

V^a{x,y)  =  I 

Y  ^'^xxy^l  “  '^xyy^2  ^'^xyy^l  “  '^yyy^2 

The  shadow  function  is  a  Morse  function  if  and  only  i/det(V^cr)  ^  0. 

5 A  Experimental  results 

The  Reeb  graph  describes  the  topological  structure  of  objects,  and  illustrates  the  topological 
changes  occured  at  singular  points  of  the  height  function  (i.e.  topological  changes  of  the  level 
sets  h^^{z)  such  as  merging  or  splitting).  The  topological  structure  offered  by  the  Reeb  graph  is 
very  useful  for  object  reconstruction  from  real  data  sets  such  as  computed  tomography  (CT)  and 
magnetic  resonance  imaging  (MRI)  usually  available  as  cross-sections.  So  we  need  to  reconstruct 
the  object  from  these  cross-sections  and  using  the  a  priori  topological  information  given  by  the 
Reeb  graph.  Furthermore,  the  Reeb  graph  describes  how  the  cells  are  glued  to  reconstruct  an 
object  surface. 

Figure  5.5  and  Figure  5.6  illustrate  the  Reeb  graph  representations  of  two  3D  models:  the  heart 
and  the  hand  models.  The  vertices  of  these  Reeb  graphs  are  the  singular  points  of  the  height 
function  in  the  light  direction  L  =  (0, 0, 1).  The  polygonal  mesh  of  the  3D  heart  model  consists  of 
861  vertices  and  1717  faces  (triangles),  whereas  the  3D  hand  object  object  is  a  laser  scanner  model 
consisting  of  38219  vertices  and  76438  faces. 

5.5  Conclusions 

Shadows  provide  perceptually  important  information  about  the  shape  of  an  illuminated  surface.  We 
have  proposed  a  new  surface  representation  function  that  provides  a  flexible  and  orientation-based 


(a)  3D  heart  model 


Figure 


(a)  3D  hand  model 


(b)  Polygonal  mesh  (c)  Reeb  graph 

5.5:  Reeb  graph  of  the  heart  model. 


Figure  5.6:  Reeb  graph  of  the  hand  model. 
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model  for  geometric  and  topological  modeling  of  surfaces.  The  geometric  and  topological  properties 
of  the  proposed  representation  function  have  been  analyzed  in  the  Morse  theory  framework. 


Chapter  6 


Geodesic  Matching  of  3D  Objects 

In  this  Chapter,  we  propose  a  new  approach  for  object  matching  based  on  a  global  geodesic  mea¬ 
sure.  The  key  idea  behind  our  methodology  is  to  represent  an  object  by  a  probabilistic  shape 
descriptor  called  geodesic  shape  distribution  [10]  that  measures  the  global  geodesic  distance  be¬ 
tween  two  arbitrary  points  on  the  surface  of  an  object.  In  contrast  to  the  Euclidean  distance  which 
is  more  suitable  for  linear  spaces,  the  geodesic  distance  has  the  advantage  to  be  able  to  capture  the 
(nonlinear)  intrinsic  geometric  structure  of  the  data.  The  matching  task  therefore  becomes  a  one¬ 
dimensional  comparison  problem  between  probability  distributions  which  is  clearly  much  simpler 
than  comparing  3D  structures.  Object  matching  can  then  be  carried  out  by  dissimilarity  measure 
calculations  between  geodesic  shape  distributions,  and  is  additionally  computationally  efficient  and 
inexpensive. 

6.1  Introduction 

Three-dimensional  objects  consist  of  geometric  and  topological  data,  and  their  compact  represen¬ 
tation  is  an  important  step  towards  a  variety  of  computer  vision  applications  including  indexing, 
retrieval,  and  matching  in  a  database  of  3D  models.  The  latter  will  be  the  focus  of  the  present 
Chapter.  There  are  two  major  steps  in  object  matching.  The  first  step  involves  finding  a  reliable 
and  efficient  shape  representation  or  descriptor,  and  the  second  step  is  the  design  of  an  appropriate 
dissimilarity  measure  for  object  comparison  between  the  shape  representations. 

Most  three-dimensional  shape  matching  techniques  proposed  in  the  literature  of  computer 
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graphics,  computer  vision  and  computer-aided  design  are  based  on  geometric  representations  which 
represent  the  features  of  an  object  in  such  a  way  that  the  shape  dissimilarity  problem  reduces  to  the 
problem  of  comparing  two  such  object  representations.  Feature-based  methods  require  that  features 
be  extracted  and  described  before  two  objects  can  be  compared.  Among  feature-based  methods, 
one  popular  approach  is  graph  matching,  where  two  objects  are  represented  by  their  graphs  com¬ 
posed  of  vertices  and  edges.  An  efficient  representation  that  captures  the  topological  properties 
of  3D  objects  is  the  Reeb  graph  descriptor  proposed  by  Shinagawa  et  al.  [67,68].  The  vertices  of 
the  Reeb  graph  are  the  singular  points  of  a  function  defined  on  the  imderlying  object  [67,68,31]. 
These  singularities  are  prominent  landmarks  and  their  detection,  recognition,  and  classification  is 
a  crucial  step  in  image  processing  and  computer  vision  [31].  Such  singularities  carry  important 
information  for  further  operations,  such  as  image  registration,  shape  analysis,  motion  estimation, 
object  recognition,  and  surface/image  evolution  [53,21,52]. 

An  alternative  to  feature-based  representations,  called  shape  distribution,  is  developed  by  Osada 
et  al  [59].  The  idea  here  is  to  represent  an  object  by  a  global  histogram  based  on  the  Euclidean 
distance  defined  on  the  siurface  of  an  object.  The  shape  matching  problem  is  then  performed  by 
computing  a  dissimilarity  measure  between  the  shape  distributions  of  two  arbitrary  objects.  This 
approach  is  computationally  stable  and  relatively  insensitive  to  noise.  Because  of  unsuitability  of 
the  Euclidean  distance  when  dealing  with  nonlinear  manifolds,  the  shape  distribution,  however, 
does  not  capture  the  nonlinear  geometric  structure  of  the  data. 

The  geodesic  shape  distribution  may  be  used  to  facilitate  representation,  indexing,  retrieval, 
find  object  matching  in  a  database  of  3D  models.  More  importantly,  the  geodesic  shape  distribution 
provides  a  new  way  to  look  at  the  object  matching  problem  by  understanding  the  intrinsic  geometry 
of  the  shape. 

Information-theoretic  measures  provide  quantitative  entropic  divergences  between  two  proba¬ 
bility  distributions.  A  common  entropic  dissimilarity  measure  is  the  Kulback-Liebler  (or  directed) 
divergence  [50]  which  has  been  successfully  used  in  many  applications  including  indexing  and  image 
retrieval  [69].  Another  entropy-based  measure  is  the  Jensen-Shannon  divergence  which  may  be  de¬ 
fined  between  an  arbitrary  number  of  probability  distributions  [51].  Due  to  this  generalization,  the 
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Jensen-Shannon  divergence  may  be  used  as  a  coherence  measure  between  any  number  of  distribu¬ 
tions  and  may  be  applied  to  a  variety  of  signal/image  processing  and  computer  vision  applications 
including  graph  matching  [39],  image  edge  detection  [37]  and  segmentation  of  DNA  sequences  into 
homogenous  domains  [64]. 

The  rest  of  this  Chapter  is  organized  as  follows.  The  next  section  is  devoted  to  the  problem 
formulation.  Section  6.3  describes  some  the  related  work  to  our  proposed  approach  for  3D  object 
matching.  In  Section  6.4,  we  describe  the  representation  step  of  our  proposed  methodology.  In 
Section  6.5,  we  present  the  Jensen-Shannon  divergence  measure  and  show  its  attractive  properties 
as  a  dissimilarity  measure  between  probability  distributions.  Section  6.6  presents  an  information 
geometric  approach  to  geodesic  shape  distributions.  In  Section  6.7,  we  provide  numerical  simula¬ 
tions  to  show  the  power  of  the  proposed  shape  measure  for  object  matching.  And  finally.  Section 
6.8  provides  some  conclusions. 

6.2  Problem  formulation 

Three-dimensional  objects  are  usually  represented  as  polygonal  or  triangle  meshes  in  computer 
graphics  and  geometric-aided  design.  A  triangle  mesh  M  is  a  pair  M  =  (V,  T),  where  V  = 
{«!, . . . ,  Vm)  is  the  set  of  vertices,  and  T  =  {Ti, . . .  ,r„}  is  the  set  of  triangles. 

In  scientific  visualization  and  analysis,  a  triangle  mesh  is  too  large  to  be  examined  without 
simplification.  One  way  to  overcome  this  limitation  is  to  represent  a  triangle  mesh  by  its  surface 
features  that  can  easily  be  computed  and  can  effectively  characterize  the  global  surface  shape.  The 
centroids  of  the  set  of  triangles  T  are  desirable  features  which  may  be  computed  efficiently  and  have 
a  global  significance  for  the  surface  shape  representation  as  illustrated  in  Figmre  6.1.  In  addition, 
there  is  a  well  defined  correspondence  between  the  centroid  and  the  region  (triangle)  firom  which 
it  is  computed  as  depicted  in  Figure  6.1.  It  is  important  to  point  out  that  centroid-based  methods 
have  been  used  in  a  variety  of  computer  vision  applications  including  clustering,  and  one  of  the 
widely  centroid-based  technique  used  for  cluster  analysis  in  the  K-mean  algorithm  [54]. 
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6.2.1  Global  shape  measure 

Let  M  =  (V,  T)  be  a  triangle  mesh.  The  centroid  Cj  of  a  triangle  Tj  is  the  mean  of  its  vertices,  that 
is,  the  point  located  at  the  center  of  the  triangle.  Note  that  the  cardinality  of  the  set  of  centroids 
C  =  {ci, . . . ,  Cn}  of  the  triangle  mesh  M  is  equal  to  the  cardinality  of  its  set  of  triangles  T . 

Unless  we  establish  a  meaningful  measure  of  distance  between  the  centroids  of  a  triangle  mesh, 
no  meaningful  exploration  of  the  underlying  structure  of  an  object  is  possible.  In  order  to  take  into 
account  the  interaction  between  the  centroids,  we  compute  a  pairwise  distance  measure  d(cf,Cj) 
from  any  centroid  Ci  to  all  the  other  centroids  Cj  G  C.  Figure  6.1  illustrates  an  arbitrary  distance 
between  two  centroids.  Notice  that  distance  d  need  not  be  an  Euclidean  metric. 


Figure  6.1:  Distance  between  two  arbitrary  centroids  of  a  3D  camel. 

To  obtain  a  global  measiure  of  the  shape  M,  we  simply  integrate  over  all  centroids.  More 
precisely,  we  define  a  function  /  :  C  C  M  R  such  that 

/(Ct)  =  1^  ^ 

where  dcj  denotes  the  area  element  that  contains  the  centroid  Cj^  that  is,  in  the  discrete  domain, 
dcj  =  area{Tj)  the  area  of  the  triangle  Tj,  and  |C|  =  (irea{Tj)  is  the  total  area  of  the  surface 
M.  The  function  /  clearly  represents  a  global  measure  or  a  distribution  of  the  shape  and  therefore 
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to  each  3D  model  M  we  will  assign  its  global  measure  / . 

The  problem  addressed  in  this  Chapter  may  now  be  concisely  described  by  the  following  state¬ 
ment:  Given  two  3D  objects  Mi  and  M2  to  be  matched,  find  their  global  measures  fi  and  /2,  and 
calculate  how  dissimilar  these  objects  are  based  on  a  predefined  dissimilarity  measure  D(/i,  /2).  In 
other  words,  the  dissimilarity  between  two  objects  measures  “how  different  they  are” ,  and  a  smaller 
value  of  D  means  that  the  two  objects  are  more  similar.  Figure  6.2  depicts  a  block-diagram  of  the 
proposed  framework. 


Figure  6.2:  Block-diagram  of  the  proposed  methodology. 
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6.2,2  Construction  of  a  measure  space 

A  measure  space  on  a  3D  model  M  is  given  by  the  triple  (C,  6,  ^),  where  B  is  a  <T-algebra  of  subsets 
of  C,  and  /r  is  a  measure  defined  on  the  set  of  triangles  T  as  /i(Tj)  =  area{Tj)  =  dcj.  Note  that 

n  n 

niC)  =  n{Tj)  =  Y  area(Tj)  <  oo, 

i=i 

hence  fi  is  a  cr-finite  measure. 

If  /i(C)  =  1,  then  fjL  is  a  probability  measure  and  therefore  we  may  define  a  random  vector 
X  :  C  such  that  X{c)  represents  the  coordinates  of  the  centroid  c  in  the  Euclidean  space. 

The  expected  value  of  X  is  given  by 

£;(X)  =  min  f  d{p,X{c)fdn  (2) 

where  d  is  a  distance  function  defined  along  the  surface  M.  This  expected  value  provides  a  nice 
statistical  interpretation  of  our  global  shape  measure  defined  in  Equation  (1). 

6.3  Related  work 

In  this  section,  we  will  review  two  representative  methods  for  object  matching  that  are  closely 
related  to  our  proposed  approach.  We  briefly  show  their  mathematical  foundations  and  edgorithmic 
methodologies  as  well  as  their  limitations. 

6.3.1  Reeb  graph 

Morse  theory  explains  the  presence  and  the  stability  of  critical  points  in  terms  of  the  topology  of 
the  underlying  smooth  manifold  [55].  Topology  is  the  property  that  determines  which  parts  of  the 
shape  of  objects  are  connected  to  which  other  parts  [38],  while  geometry  determines  where,  in  a 
given  coordinate  system,  each  part  is  located  [29].  The  basic  principle  is  that  the  topology  of  a 
manifold  is  very  closely  related  to  the  critical  points  of  a  smooth  function  on  that  manifold. 

An  interesting  concept  related  to  Morse  theory  and  very  useful  to  analyze  a  surface  topology 
is  the  Reeb  graph.  The  latter  is  defined  as  a  quotient  space  M/~  with  an  equivalence  relation 
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given  by  p  ~  g  if  and  only  if  h{p)  —  h{q)  and  p,  q  belong  to  the  same  connected  component  of 
/i“^(/i(p)),  where  /i  :  M  — >  E  is  the  height  function  such  that  h{x^y,z)  =  z  for  all  (x^y^z)  €  M. 
In  particular,  each  connected  component  is  represented  by  a  point  in  the  Reeb  graph  as  illustrated 
in  Figure  6.3.  The  left  figure  shows  a  torus  with  the  critical  points  of  its  height  function  (Morse 
function).  The  figure  in  the  middle  illustrates  the  geometric  featmres  of  the  torus  represented  by 
cross-sections  along  its  height.  The  right  figure  shows  the  topological  features  represented  by  the 
Reeb  graph.  By  taking  an  appropriate  number  of  cross-sections  and  smooth  interpolation  between, 
Shinagawa  et  al  [67, 68]  proposed  a  Reeb  graph  based  approach  or  so-called  homotopy  model  for 
object  reconstruction. 


Figure  6.3:  Reeb  graph  representation  of  a  torus. 

Figure  6.4  illustrates  the  Reeb  graph  representation  of  a  3D  hand  model.  The  vertices  of  the 
Reeb  graph  are  the  singular  points  of  the  height  function.  The  polygonal  mesh  of  the  3D  hand 
object  object  is  a  laser  scanner  model  consisting  of  38219  vertices  and  76438  faces. 

Reeb  graph  has  a  nice  mathematical  definition  that  makes  it  very  attractive  from  a  theoretical 
point  of  view.  This  representation,  however,  is  not  rotationally  invariant.  This  limitation  lead 
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(a)  3D  hand  model 


Figure  6.4:  Reeb  graph  of  the  hand  model. 


Hilaga  et  al  to  develop  a  geodesic-based  Reeb  graph  technique  [40].  In  this  approach  a  multires¬ 
olution  Reeb  graph  is  computed  efficiently  and  a  similarity  distance  is  calculated  to  compare  two 
Reeb  graphs. 


6.3.2  Shape  distribution 

Recently,  Osada  et  al  [59]  proposed  a  shape  distribution  based  approach  for  three-dimensional  object 
matching.  The  key  idea  is  to  compute  the  Euclidean  distance  between  all  pairs  of  random  points 
on  the  surface  to  obtain  the  so-called  D2  shape  histogram.  Given  a  triangle  Tj  =  {vji,Vj2,Vj2}, 
each  r2indom  point  is  generated  as 


Pj  =  (1  -  \/ri)vjx  -b  -  r{)Vj2  +  v/nravp, 

where  ri  and  r2  are  pseudo-random  numbers  between  zero  and  the  total  cumulative  area.  Then,  the 
comparison  of  objects  is  carried  out  by  computing  a  dissimilarity  measure  between  their  D2  shape 
distributions.  Figure  6.5  illustrates  an  ellipsoid  emd  its  D2  shape  distribution.  The  main  drawback 
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of  the  shape  distribution  approach,  which  is  based  on  the  Euclidean  distance,  is  its  inability  to 
capture  the  nonlinear  structure  of  the  data. 


(a)  Ellipsoid  (b)  D2  shape  distribution 

Figure  6.5:  D2  shape  distribution  of  an  ellipsoid. 

6.4  Proposed  methodology 

The  goal  of  our  proposed  approach  may  be  described  as  follows:  Given  two  3D  objects  Mi  and 
M2  to  be  matched,  find  their  global  measures  or  shape  descriptors  fi  and  /2,  and  calculate  how 
dissimilar  these  objects  using  a  dissimilarity  measure  D(fi,  that  has  to  be  quantified.  The  basic 
idea  behind  the  shape  descriptor  is  to  characterize  a  3D  object  with  a  one-dimensional  function 
which  will  help  us  discriminate  between  objects  in  a  database  of  3D  models. 

6.4.1  Global  geodesic  shape  function 

The  Reeb  graph  concept  has  been  shown  to  be  very  effective  in  modeling  3D  objects  based  on 
cross  sections  such  as  MRI  or  CT  images.  It  is  more  appropriate  to  modeling  applications  where 
the  height  is  of  special  interest  such  as  terrain  imaging.  The  height  function,  however,  has  some 
limitations  as  an  object  signature  for  matching,  indexing,  or  retrieval  of  arbitrary  3D  objects.  The 
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main  reason  is  that  the  height  function  is  not  rotationally  invariant.  To  overcome  these  limitations, 
we  propose  a  global  geodesic  function  defined  on  the  object  surface  as  follows.  Let  Ci  and  Cj  be 
two  points  (centroids)  on  a  manifold  M.  The  geodesic  distance  g{ci^Cj)  between  Ci  and  Cj  is  the 
shortest  length  L{'y)  =  ||7'(t)||d<  of  a  smooth  curve  7  :  [a,6]  — ►  M  such  that  7(a)  =  Ci  and 
7(6)  =  Cj,  The  geodesic  distance  may  be  locally  viewed  as  the  Euclidean  d^(ci,Cj)  =  ||ct  --  Cj||, 
and  is  hence  clearly  invariant  to  rotation  and  translation. 

Inspired  by  the  geodesic-based  representation  for  3D  topology  matching  proposed  by  Hilaga  et 
aL  [40] ,  we  define  a  global  shape  function  f  :  C  expressed  in  terms  of  a  rotationally  invariant 
(square)  geodesic  distance  as  follows 

/(ci)  =  ^  (3) 

The  primary  motivation  behind  the  geodesic  distance  is  of  overcoming  the  limitations  of  the  Eu¬ 
clidean  distance  which  by  virtue  of  its  linearity  in  nature  cannot  account  for  nonlinear  structures 
in  the  underlying  object. 

Unlike  the  Euclidean  distance  which  is  basically  a  straight  line  between  two  points  in  3D  space, 
the  geodesic  distance  captures  the  global  nonlinear  structure  and  the  intrinsic  geometry  of  the  data 
as  illustrated  in  Figure  6.6.  This  clearly  shows  that  the  Euclidean  distance  between  two  arbitrary 
points  in  a  nonlinear  manifold  is  just  a  straight  segment  connecting  two  points  and  does  not  reflect 
the  nonlinear  structure  of  the  object,  whereas  the  geodesic  distance  which  is  the  shortest  curve 
along  the  manifold  connecting  both  points  clearly  captures  the  intrinsic  geometry  of  the  object. 

Geodesic  distance  calculation 

Given  a  set  of  centroids  C  —  {ci, . . . ,Cn}  of  a  3D  object  represented  by  a  triangle  mesh  M,  the 
geodesic  distance  calculation  is  based  on  a  similar  approach  used  for  computing  the  isometric  feature 
mapping  (Isomap)  for  multidimensional  scaling  on  nonlinear  manifolds  [71].  The  algorithm  has  two 
main  steps: 

(i)  Construct  a  neighborhood  graph  by  connecting  a  given  centroid  to  its  fc-nearest  neigh¬ 
bors,  and  link  these  neighboring  centroids  by  edges  with  weights  equal  to  the  Euclidean 
distances. 
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Figure  6.6:  Euclidean  vs.  geodesic  distance  on  a  nonlinear  manifold. 


(ii)  Compute  the  geodesic  distances  (shortest  paths)  between  all  pairs  of  n  points  in  the 
constructed  graph  using  Dijkstra’s  or  Floyd’s  algorithm 
Note  that  Step  (i)  may  be  alleviated  by  choosing  a  random  subset  of  C  in  order  to  reduce  the 
computational  complexity  of  the  geodesic  calculation. 

Prom  Equation  (3),  it  is  clear  that  a  discrete  form  of  the  geodesic  shape  function  can  be  written 


f{Ci)  = 


{Ga)i  iGa)i 


where  G  =  (gfj)  is  the  (square)  geodesic  distance  matrix  of  size  n  x  n,  and  a  = 
is  an  n  X  1  vector  of  triangle  areas,  i.e.  aj  =  area{Tj),  and  \C\  =  ~  ll®lli  total 

area.  The  geodesic  distance  matrix  G  =  (5?  )  is  symmetric  with  zeros  in  the  diagonal,  and  positive 
off-diagonal  elements 
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TViangle  area  calculation 


Denote  by  {v\,V2,vz}  the  vertices  of  an  arbitrary  triangle  T  of  a  given  triangle  mesh  M.  Using 
Newell  method,  the  area  of  the  triangle  T  can  be  computed  as 


a  =  area{T)  = 


m 

2 


where  N  =  (iVi,  iV2, ,  Ns)  is  the  triangle  normal  vector  with  coordinates  given  by 

d 

^  ynext(t))('^i  “I"  ^next(i)) 

i=l 
d 

N2  =  ^  ~  '^next(i))(^i  ^next(i)) 

t=l 

d 

jV3  =  ^  ^next(t))(yt  “1“  ynext(t)) 

t=l 

and  {xi^  2/i,  Zi)  are  the  coordinates  of  eaeh  vectex  Vi  (with  dimension  d  =  3)  of  a  triangle  T.  Note 
that  rnext(t)  =  (a^next(i),ynext(f),^next(t))  denotes  the  next  vertex  in  the  list  after  vu  taking  into 
account  that  vi  follows  the  last  vertex  Vd- 


6.4.2  Global  geodesic  shape  distribution 

Note  that  the  geodesic  shape  function  can  be  expressed  as  a  geodesic  shape  vector  X  =  {Xi, . . . ,  Xn}, 
where  Xi  =  /(cj).  This  vector  may  be  viewed  as  a  shape  descriptor  that  may  be  used  for  3D  shape 
comparison. 

Assume  that  the  geodesic  shape  vector  X  of  an  object  M  is  a  random  sample  with  a  common 
(unknown)  probability  density  function  p.  A  common  approach  to  approximate  the  probability 
density  function  p  is  through  the  kernel  density  estimation  which  is  an  important  data  analytic 
tool  that  provides  a  very  effective  way  of  showing  structure  in  a  data  set  [73] .  The  kernel  density 
estimator  p  is  given  by 

p) 

where  K  is  the  Gaussian  kernel,  and  h  is  the  bandwidth  or  window  width  to  be  estimated.  A  good 
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Figure  6.7:  Effect  of  the  bandwidth  parameter  h. 

selection  rule  of  this  bandwidth  is  given  by 

^  [  24S  R{K)  1^/®  . 

^“[35/i2WnJ 

where  RiK)  =  J K{t)'^dt,  P2{K)  =  /  t'^K{t)dt,  and  d  =  medj{|Xj  -  medi{Xi}|}  is  the  median 
absolute  deviation.  The  effect  of  the  bandwidth  parameter  h  is  illustrated  in  Figure  6.7. 

Hence  to  each  3D  object  represented  by  a  triangle  mesh  M,  we  associate  a  kernel  density  p 
which  we  will  refer  to  as  a  geodesic  shape  distribution,  and  it  is  computed  using  the  algorithm 
depicted  in  Figure  6.8.  This  probabilistic  shape  descriptor  represents  an  object  information  and 
will  be  used  in  our  matching  experiments.  Figure  6.9  depicts  a  3D  model  of  a  tank  and  its  geodesic 
shape  distribution. 

6.4.3  Properties  of  geodesic  shape  signature 

In  addition  to  its  rotational,  translational  and  scale  invariance,  the  geodesic  shape  signature  is 
also  robust  to  resampling  and  simplification  as  illustrated  in  Figures  6.10  and  6.11.  Note  that 
for  triangulation,  we  use  the  barycentric  subdivision  shown  in  the  top  row  of  Figure  6.10.  This 
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subdivision  technique  consists  of  introducing  a  new  vertex  at  the  center  of  each  triangle  and  a  new 
vertex  at  the  midpoint  of  each  edge  and  drawing  edges  from  the  centroid  of  the  triangle  to  each  of 
the  new  midpoint  vertices  and  to  the  original  vertices. 


Figure  6,10:  Robustness  and  invariance. 


In  order  to  compare  two  geodesic  shape  distributions  and  hence  to  measure  the  performance  of 
the  proposed  scheme,  we  will  describe  in  the  next  section  an  information-theoretic  distance  that 
quantifies  the  difference  between  two  3D  shapes  through  their  probabilistic  shape  descriptors. 


6.5  Probabilistic  dissimilarity 

Let  Ml  and  M2  be  two  3D  objects  with  geodesic  shape  distributions  p  and  q  respectively.  In¬ 
formation  theoretic  measures  provide  quantitative  entropic  divergences  between  two  probability 
distributions.  A  common  entopic  dissimilarity  measure  is  Kulback-Liebler  (KL)  divergence  /C  de¬ 
fined  as 

ACtf,«  =  /p(x)log,Mir  =  E{log^}, 
where  E{-}  denotes  the  expected  value  with  respect  to  p(x). 
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(c)  Ms  (d)  geodesic  shape  distributions 

Figure  6.11:  Robustness  and  invariance  (cont.). 


The  KL  dissimilarity  measure,  however,  is  non-symmetric,  unbounded,  and  undefined  if  p  is 
not  absolutely  continuous  with  respect  to  g  [37].  To  overcome  these  limitations,  we  use  the  Jensen- 
Shannon  divergence  D  given  by 

DM)  = 

+  H(p)  +  Hiq)^ 

where  H(p)  =  —  f  p(x)  log2  p(x}  dx  is  the  differential  entropy,  which  corresponds  to  Shannon’s 
entropy  in  the  discrete  domain.  Shannon’s  entropy  is  a  measure  of  uncertainty,  dispersion,  in¬ 
formation,  and  randomness.  The  maximum  uncertainty  or  equiveJently  minimum  information  is 
achieved  by  the  uniform  distribution.  Hence,  we  can  think  of  the  entropy  as  a  measure  of  uniformity 
of  a  probability  distribution.  Consequently,  when  uncertainty  is  higher  it  becomes  more  difficult  to 
predict  the  outcome  of  a  draw  from  a  probability  distribution. 
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The  Jensen-Shannon  divergence  is  a  statistical  distance  that  is  very  useful  in  quantifying  dif¬ 
ferences  between  probability  distributions  or  densities.  In  other  words,  this  dissimilarity  measure 
quantifies  differences  in  shape  between  two  arbitrary  objects.  Unlike  the  Kullback-Leibler  diver¬ 
gence,  the  Jensen-Renyi  divergence  has  the  advantage  of  being  symmetric,  always  defined,  and 
generalizable  to  any  arbitrary  number  of  probability  distributions,  with  a  possibility  of  assigning 
weights  to  these  distributions  [51].  Figure  6.12  shows  a  three-dimensional  graph  and  a  contour  plot 
of  the  Jensen-Shannon  divergence  between  two  discrete  Bernoulli  distributions. 


Figure  6.12:  (a)  3D  plot  and  (b)  contour  plot  of  the  Jensen-Shannon  divergence. 

The  following  result  establishes  the  convexity  of  the  Jensen-Shannon  divergence. 

Proposition  6.1  The  Jensen-Shannon  divergence  D{p^q)  is  a  convex  function  of  p  and  q. 

In  addition  to  its  convexity  property,  the  Jensen-Shannon  divergence  is  shown  to  be  an  adapted 
measure  of  disparity  among  probability  distributions.  Using  the  theory  of  majorization,  it  can  be 
shown  that  the  Jensen-Shannon  divergence  is  bounded,  and  its  upper  bound  is  achievable. 

Proposition  6.2  The  Jensen-Shannon  divergence  between  two  geodesic  shape  distributions  p  and 
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q  is  upper  bounded 


D(p,q)  <  log2(2)  =  1. 


6.6  Information-geometric  approach  to  geodesic  shape  distribu¬ 
tions 


The  Jensen-Shannon  divergence 


Dip,q)  =  2 


K 


satisfies  the  triangle  inequality  property 

DipuP2)  +  Dip2,P3)  >  Dipups) 


if  and  only  if 

H  (51^)  +  H  j 

While  the  Jensen-Shannon  divergence  is  not  a  metric,  it  can  be  shown  that  its  square  root 
\/D{',  •)  is  a  metric  between  probability  distributions. 


6.6.1  Statistical  manifolds 

A  Riemannian  manifold  is  a  differentiable  manifold  equipped  with  a  positive  definite  inner  product 
<  •,  •  >x'-  TxM.  X  TxM.  — >  K.  The  collection  of  all  these  inner  products  is  called  the  Riemannian 
metric.  An  example  of  such  a  metric  is  the  first  fundamental  form  derived  in  the  Appendix. 

Statistical  Manifolds  are  differentiable  meinifolds  such  that  each  point  can  be  identified  with  a 
probability  density  with  respect  to  a  given  measure,  and  a  family  of  distributions  correspond  to 
a  set  of  points  which  form  a  manifold.  The  theory  of  statistical  manifolds  also  called  information 
geometry  was  originally  proposed  by  Rao  [63]  who  considered  a  parametrized  statistical  model  as  a 
Riemannian  manifold  with  the  metric  tensor  given  by  the  Fisher  information  metric.  This  metric 
defines  a  notion  of  a  distance  between  two  probability  distributions.  In  other  words,  it  measures 
“how  far  apart  are  these  distributions  ?”.  A  good  reference  to  information  geometry  is  the  book 
by  Amari  [19]  who  introduced  a-connections  and  showed  how  they  relate  to  asymptotic  inference. 
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Let  M  be  a  differentiable  manifold  representing  a  statistical  model  M  =  {p(x;  0)  :  ff  €  0} 
of  probability  distributions  p(x;  0)  pareunetrized  by  a  real-valued  vector  0.  In  other  words,  Af  is 
a  parametric  surface  which  can  be  represented  as  a  Monge  patch  defined  by  r  :  0  — >  Af  such 
that  r(0)  =  {0,p{x]0)).  Note  that  the  patch  r  covers  all  M,  that  is,  r(0)  =  Af.  It  is  worth 
pointing  out  that  the  Monge  patch  r  :  0  — +  Af  of  a  smooth  probability  density  p(x;  •)  :  0  — >  K  is  a 
diffeomorphism  because  it  is  a  smooth  bijection,  and  its  inverse  r~^  is  the  restriction  to  Af  of  the 
smooth  projection  tt  :  0  x  K  0,  that  is  r~^  =  i:\m' 

The  Kullback-Leibler  divergence  two  points  p{x\  0)  and  p(x;  0')  in  Af  is  given  by 


AC«-, «')  =  /  p(., «)  log  ^  <Jx  =  E  {log  ^  , 

where  E{-}  denotes  the  expected  value  with  respect  to  p{x;  0). 

When  0  and  0'  are  infinitesimally  close  to  each  other  (i.e.  0*  =  0  +  e  with  e  sufficiently  small), 
it  can  be  shown  that 


IC{0,0')  =  1(0'-  0fT.{0W  -  ^)  +  Oi\\0'  - 

A 


where  S(0)  =  {crij{0))  is  the  Fisher  information  metric  tensor  given  by 

(rtj{0)  =  j p{x\  0)di  logp(x;  0)dj  logp(x;  0)  dx  =  E{di  logp(x;  0)dj  logp(x;  0)}, 
and  di  logp(s;  0)  denotes  the  partial  derivative  with  respect  to  the  i-th  component  of  the  vector  0. 


6.6.2  Geodesic  shape  manifold 


Let  M  =  (V, be  a  triangle  mesh,  where  V  =  {ui,---,v,n}  is  the  set  of  vertices,  and  T  = 
is  the  set  of  triangles.  Denote  by  a  the  total  area  of  the  triangle  mesh,  that  is 

«  =  E”=i  «»’ea(Tj). 

To  apply  information  geometry  to  our  proposed  geodesic  shape  distribution  and  using  Equation 
(4),  we  may  rewrite  the  geodesic  kernel  density  in  parametric  form  as 


p{x\  0) 


»)), 
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where  0  =  (n,  a)  defines  the  distribution  parameters:  the  number  of  triangles  n  and  the  surface 
total  area  a.  The  3D  object  matching  by  the  information  geometric  approach  may  be  stated 
as  follows:  Given  a  database  of  3D  models  {A^i, the  first  step  is  to  compute  the 
corresponding  parametrized  probability  shape  distributions  {p{x\  6i),p(rr;  62)5  •  •  •  fhat  is 

each  3D  model  can  be  viewed  as  a  point  in  a  set  <S  =  {p{x\0)  :  6  €  N  x  E4.}  embedded  in  the 
3D  Euclidean  space  as  displayed  in  Figure  6.13.  This  set  S  of  geodesic  shape  distributions  forms  a 
statistical  model  that  carries  the  structure  of  a  smooth  manifold,  and  will  be  referred  to  as  geodesic 
shape  manifold.  Note  that  the  parameter  6  =  (n,  a)  plays  the  role  of  coordinates  of  a  geodesic 
shape  density  p{x]  6)  E  S, 


Figure  6.13:  Illustration  of  geodesic  shape  statistical  manifold. 

6.7  Experimental  results 

Object  matching  experiments  were  performed  using  a  database  of  3D  models  collected  online. 
Each  model  is  represented  as  a  triangle  mesh.  We  conducted  four  sets  of  experiments.  The  first  set 
consists  of  3D  airplane  models  as  shown  in  Figure  6.14,  and  the  second  set  consists  of  3D  tanks  as 
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illustrated  in  Figure  6.15. 


Figure  6.14:  First  set  of  experiments:  3D  airplanes. 

The  third  set  deals  with  objects  that  are  topologically  equivalent  to  a  sphere  (i.e.  with  genus 
equal  to  zero)  as  shown  in  Figure  6.16.  The  numerical  results  using  the  Jensen-Shanon  dissimilarity 
measure  are  depicted  in  Table  6.1  where  the  grayscale  colorbar  displays  the  grayscale  colormap  of 
this  dissimilarity  matrix.  This  grayscale  colormap  ranges  from  white  (maximum  similarity)  to  black 
(maximum  dissimilarity),  and  passes  through  the  gray  colors  indicating  the  values  of  the  matching 
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Figure  6.15:  Second  set  of  experiments:  3D  tanks. 

algorithm.  Note  that  the  minimum  dissimilarity  rate  is  about  9%,  that  is  the  matching  rate  is 
about  81%. 

The  fourth  set  of  experiments  is  similar  to  the  third,  except  that  the  underlying  objects  are 
topologically  different  from  than  the  ones  in  the  third  set  of  the  experiments.  Figure  6.17  shows 
a  set  of  objects  with  genus  equal  to  one.  Matching  is  achieved  by  the  minimum  Jensen-Shannon 
distance  computations  as  illustrated  in  Table  6.2.  Note  that  the  minimum  dissimilarity  rate  is 
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about  2%,  that  is  the  matching  rate  is  about  98%. 
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6.8  Conclusions 

In  this  Chapter,  we  proposed  an  new  methodology  for  object  matching.  The  key  idea  is  to  encode 
a  3D  shape  into  a  ID  geodesic  shape  distribution.  Object  matching  is  then  achieved  by  calculat¬ 
ing  an  information-theoretic  measure  of  dissimilarity  between  the  probability  distributions.  That 
is,  the  dissimilarity  computations  are  carried  out  in  a  low-dimensional  space  of  geodesic  shape 
distributions.  The  main  advantages  of  the  proposed  approach  are: 

•  The  geodesic  distance  captures  the  intrinsic  geometry  of  the  data 

•  The  approach  is  simple  and  computationally  inexpensive 

•  The  simulations  results  indicate  the  suitability  of  the  proposed  technique  for  object 
matching 

For  future  work,  it  would  be  of  interest  to  incorporate  topology  into  the  proposed  methodology 
through  Morse  singularities  of  the  global  geodesic  shape  function.  Finally  we  note  that  while  the 
experimental  results  presented  in  this  section  are  very  promising,  significant  additional  performance 
gains  are  still  possible.  For  example,  our  current  way  of  selecting  centroids  as  landmarks  is  rather 
one  of  many  possible  options  and  by  no  means  the  best  option,  and  a  multiresolution  geodesic 
shape  distribution  may  also  provide  better  key  to  landmarks. 


Chapter  7 


Distance  Function-based  Object  Recognition 

In  this  Chapter  we  propose  a  distance  function-based  approach  to  topological  modeling  of  3D 
objects.  Despite  the  theoretical  nature  of  the  results  presented  in  this  Chapter,  the  key  idea  is  to 
identify  and  encode  regions  of  topological  interest  of  a  3D  object  in  the  Morse-theoretic  framework. 
The  main  motivation  behind  using  the  distance  function  is  its  rotational  invariance  which  makes 
it  more  adapted  to  object  recognition  than  the  Morse  height  function.  We  prove  that  a  surface 
may  be  reconstructed  from  its  intersections  with  concentric  spheres  centered  at  the  barycenter  of 
the  underlying  surface.  The  topological  changes  in  the  surface  occur  as  we  increase  the  value  of 
the  sphere  radius.  At  singular  points,  the  level  curves  of  the  distance  function  may  split  or  merge 
which  indicate  topological  changes.  We  also  show  that  when  a  surface  is  embedded  in  a  sphere,  the 
height  function  and  the  distance  function  are  equivalent  in  a  Morse-theoretic  setting,  that  is  both 
functions  have  the  same  singularities. 

7.1  Introduction 

In  computer  graphics  applications,  one  is  typically  interested  in  locating  geometric  regions  of  topo¬ 
logical  nature  on  a  surface.  The  simplest  non-trivial  regions  are  areas  with  genus  equal  to  one. 
Such  regions  are  handles.  As  mentioned  in  Chapter  1,  a  handle  intuitively  corresponds  to  its  def¬ 
inition  per  se.  For  example,  a  coffee  mug  with  a  handle  has  genus  equal  to  one.  Our  challenge  in 
this  research  effort  is  to  present  computational  topology  algorithms  which  are  adapted  to  discrete 
surfaces  and  which  simultaneously  account  for  the  geometry  of  a  surface.  The  main  objectives  of 
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this  chapter  consist  of: 

•  identifying  regions  of  topological  interest,  that  is  where  topology  changes  on  a  surface 
occur.  Identifying  the  topology  of  a  surface  is  tightly  related  to  Morse  theory  which 
establishes  a  relationship  of  critical  points  of  a  smooth  function  defined  on  a  smooth 
manifold  to  its  connectivity. 

•  coding  a  surface  topology  into  a  Reeb  graph.  The  nodes  of  this  graph  represent  the 
critical  points  of  a  function  defined  on  the  surface,  and  the  edges  in  the  graph  represent 
the  connected  components  of  the  surface  between  critical  points. 

7.2  Topology  identification 

Topology  is  a  branch  of  mathematics  dealing  with  qualitative  questions  about  geometrical  struc¬ 
tures.  We  do  not  ask:  how  big  is  it?  but  rather:  does  it  have  any  holes  in  it?  Is  it  all  connected 
together,  or  can  it  be  separated  into  parts?  Geometry,  on  the  other  hand,  deals  with  measuring  and 
computing  lengths,  areas,  volumes,  angles  etc.,  and  that  is  actually  where  the  word  “geo-metry” 
comes  from.  The  subject  of  topology  is  concerned  with  those  features  of  geometry  which  remain 
imchanged  after  twisting,  stretching  or  other  deformations  of  a  geometrical  space.  It  includes  such 
problems  as  distinguishing  knots  and  classifying  surfaces.  One  of  the  key  tools  used  to  study  the 
topology  of  spaces  is  Morse  theory  which  is  the  study  of  the  relationship  between  functions  on 
a  space  and  the  shape  of  the  space.  Although  Morse  theory  can  be  applied  to  spaces  of  infinite 
dimension,  we  are  particularly  interested  in  the  application  of  Morse  theory  to  2-manifolds.  Based 
on  the  calculus  of  variations,  Morse  theory  draws  a  relationship  between  critical  points  of  a  smooth 
function  defined  on  a  smooth  manifold  and  the  global  topology  of  that  manifold.  To  better  under¬ 
stand  Morse  theory,  we  start  by  briefly  introducing  some  basic  definitions  in  differential  geometry 
and  topology. 

Definition  7.1  An  abstract  surface  or  2-manifold  is  a  topological  space  M  such  that  each  point  of 
M  has  a  neighborhood  U  in  M  homeomorphic  to  an  open  2-disk  in  R^. 
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In  other  words,  a  2-manifold  is  locally  homeomorphic  to  an  open  disk  Homeomorphism  is  a 
continuous  function  defined  between  two  spaces  which  is  bijective  and  also  has  a  continuous  inverse. 

If  M  is  a  2-manifold,  then  we  can  find  a  countable  system  of  open  sets  Ui  and  homeomorphisms 
cpi  :  Ui  such  that  M  =  UiC/i.  This  homeomorphism  is  illustrated  in  Figure  7.1.  A  collection 

of  charts  is  called  an  atlas,  A  2-manifold  M  may  be  embedded  in  meaning  that  it  has  no 
self-intersections. 


Figure  7*1:  Definition  of  a  2-manifold. 

A  2-manifold  is  a  surface  where  the  local  area  around  every  point  on  the  surface  is  Euclidean, 
meaning,  around  each  point  the  surface  appears  to  be  nearly  fiat.  The  world  around  us  is  an 
excellent  example  of  a  2-manifold.  Manifolds  are  a  preferable  surface  representation  because  the 
surface  can  be  divided  into  charts  which  allow  2-manifolds  embedded  in  3D  to  be  fiattened  into 
a  two  dimensional  domain  (through  parametrization).  Surfaces  used  in  computer  graphics  are 
typically  oriented,  this  refers  to  the  fact  that  the  surface  has  two  sides.  For  example,  a  sphere  has 
two  sides,  while  a  Mobius  strip  has  only  one  side.  Another  attribute  of  surfaces  is  whether  the 
surface  is  closed  or  with  boundary.  This  refers  to  the  number  of  open  boundary  components  of  a 
surface.  For  example,  an  egg  shell  is  closed  but  once  it  has  been  cracked  open,  it  becomes  a  surface 
with  boundaries. 
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We  refer  to  a  surface  as  a  smooth  and  compact  2-manifold  without  boundary  and  possibly 
embedded  in  the  Euclidean  space  Mathematically,  surfaces  are  often  conceived  of  as  continuous 
and  smooth,  i.e.,  one  that  has  a  sufficient  number  of  partial  derivatives.  Smooth  often  refers 
to  a  surface  with  infinitely  many  partial  derivatives,  but  in  practice  second-order  derivatives  are 
sufficient.  In  computer  graphics  we  operate  in  a  discrete  setting,  where  only  a  finite  number  of 
samples  are  used  to  represent  a  surface.  These  surfaces  are  often  continuous,  but  are  only  piece- 
wise  linear  and  are  represented  only  by  a  discrete  set  of  points  which  are  connected  together  as 
triangles  or  polygons. 

7.2.1  Singular  points 

Let  :  M  — >  R  be  a  real- valued  function  defined  on  a  smooth  manifold  M  C  R^.  The  function  (p 
is  smooth  if  the  composition  function  tpo  p, :  U  is  smooth  (in  the  ordinary  Euclidean  sense) , 
where  p  is  &  smooth  regular  parametrization  of  M  (i.e.  p  :  U  G  R^  — »  R^).  A  point  Pq  on  M  is  a 
singularity  or  critical  point  of  p  iipo  =  p{xo,  yo),  for  some  (so,  J/o)  €  U ,  and  the  gradient  oi  cpo  p 
at  {xo,yo)  vanishes,  i.e.  V{(po  p(xo,yo))  =  0. 

A  singularity  Pq  is  nondegenerate  if  the  Hessian  matrix  V^(v7  o  p{xo,  yo))  is  nonsingular.  Oth¬ 
erwise  this  singularity  Pq  is  said  to  be  degenerate. 

7.2.2  Morse  function 

Morse  theory  explains  the  presence  and  the  stability  of  critical  points  in  terms  of  the  topology  of 
an  underlying  smooth  manifold.  Topology  is  the  property  that  determines  which  parts  of  an  object 
are  connected  to  which  other  parts  [38],  while  geometry  determines  where,  in  a  given  coordinate 
system,  each  part  is  located  [29].  The  basic  principle  is  that  the  topology  of  a  manifold  is  very 
closely  related  to  the  critical  points  of  a  smooth  function  on  that  manifold. 

Morse  proved  a  major  result  which  generalizes  the  straightforward  result  that  the  lowest-order 
nonvanishing  term  in  a  Taylor  series  describes  the  local  behavior  of  a  smooth  function  of  a  single 
variable  to  functions  of  many  variables. 

Definition  7.2  A  smooth  function  ^  :  M  — >  R  on  a  smooth  manifold  M  is  called  a  Morse  function 
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if  all  its  singularities  are  nondegenerate. 

Examples  of  nondegenerate  singular  points  are  shown  in  Figure  7.2. 


Figure  7.2:  Critical  points. 

Nondegenerate  singularities  are  isolated,  that  is,  there  ceinnot  be  a  sequence  of  nondegenerate 
singularities  converging  to  a  nondegenerate  singularity  p.  In  other  words,  there  is  no  other  point 
in  the  neighborhood  of  p  that  is  singular.  This  fact  follows  from  the  following  Morse’s  lemma. 

Lemma  7.3  Ifip  :  M  — ►  R  has  a  nondegenerate  singularity  atpQ  €  M,  then  there  exists  (xq,  yo)  €  f2 
such  that  Pq  =  r{xo,yo),  and  (p  has  the  following  representation 

(p(p)  =  <po  r(xo,  yo)  ±x^±y^  =  (p(po)  ±x^± 

for  all  p  =  r(x,  y)  e  M,  where  r  is  a  regular  smooth  path. 

Note  that  the  only  nondegenerate  singularities  are  the  minimum,  maximum  and  saddle  points.  By 
decomposing  a  smooth  manifold  along  these  singularities,  its  global  shape  and  topology  is  revealed. 
Morse  theory  also  presents  methods  to  classify  critical  points.  Specifically,  by  examining  the  number 
of  negative  eigenvalues  of  the  Hessian,  the  critical  point  can  be  indexed.  That  is,  a  minimum  has 
zero  negative  eigenvalues,  a  saddle  point  has  one,  and  a  maximum  has  two  negative  eigenvalues. 
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This  analysis  corresponds  to  the  fact  that  a  minimum  has  no  downhill  sides,  while  an  isolated 
saddle  point  has  two  downhill  sides,  one  parallel  to  the  direction  of  the  eigenvector  associated  with 
the  negative  eigenvalue  and  one  anti-parallel.  A  maximum  has  downhill  sides  associated  with  both 
directions  of  both  eigenvectors. 

7.2.3  Sard’s  theorem 

Let  /  :  M  — >  M  be  a  smooth  function.  A  point  p  is  called  a  regular  point  of  /  if  the  differential 
df  :  TpM.  ^  R  is  surjective,  that  is,  the  Jacobian  matrix  (3  x  1  in  this  case)  has  rank  equal  to 
dim(R)  =  1.  Otherwise,  the  point  p  is  called  a  critical  point.  Denote  by  Crit(/)  the  set  of  critical 
points  of  /. 

Theorem  7.4  (Sard)  The  set  /(Crit(/))  of  critical  values  of  f  has  measure  zero  in  M.  (in  the 
sense  of  Lebesgue  measure). 

Corollary  7.5  The  set  R  —  /(Crit(/))  of  regular  values  of  f  is  dense  in  R. 

Note:  /  can  be  defined  between  two  smooth  memifolds  with  arbitrary  dimensions,  i.e.,  /  :  M  — M. 

Definition  7.6  A  smooth  map  /  :  M  — >  M  is  called  an  immersion  if  at  any  point  p  €  M,  the 
differential  df  :  2pM  — >  r^(p)M  is  injective,  i.e.,  no  nonzero  vector  maps  to  zero.  If,  moreover,  f 
is  a  homeomorphism  when  considered  as  a  map  from  M  to  /(M),  we  say  that  f  is  an  embedding. 

7.2.4  Height  function 

A  classic  result  from  Morse  theory  is  that  given  a  closed  surface  M  and  a  Morse  function  /  :  M  — >  R, 
we  can  show  that  if  this  function  has  only  two  non-degenerate  critical  points  then  M  is  topologically 
equivalent  to  a  sphere.  For  example,  a  tjq)ical  Morse  function  is  a  height  function,  and  if  we  consider 
such  a  function  defined  on  a  sphere,  we  can  identify  two  critical  points  which  corresponding  to  the 
mayimiiTTi  and  minimum  at  the  north  and  south  pole  of  the  sphere  (see  Figure  7.2).  More  precisely, 
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the  the  Euler  characteristic  of  a  surface  is  defined  in  terms  of  the  number  of  critical  points  as 
follows: 

^  z=  ^minima  —  ^saddlepoints  +  ifmaxima. 

Another  essential  result  from  Morse  theory  shows  that  between  critical  points  the  topology  of  the 
manifold  is  guaranteed  not  to  change  (called  the  deformation  lemma). 

To  further  illustrate  the  relationship  of  critical  points  and  the  global  topology  of  a  smface, 
consider  the  following  geometric  interpretation.  Given  a  Morse  function,  /  :  M  — >  R  which  is  a 
height  function  which  may,  for  example,  define  parallel  planes.  A  height  function  in  the  z-direction 
on  a  smooth  manifold  M  is  a  real- valued  function  /i  :  M  R  such  that  h{x,y,z)  =  z  for  all 
(x,  y,  z)  €  M.  Hence,  h  is  the  orthogonal  projection  with  respect  to  the  z-axis.  Figure  7.3  shows 
a  2D  manifold  (a  double  torus)  and  the  corresponding  critical  points  of  its  height  function.  These 
singular  points  are  all  nondegenerate. 
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Figure  7.3:  A  3-D  object  and  the  critical  points  of  its  height  function. 
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Now  imagine  a  torus  standing  on  its  end  (see  Figure  7.2).  When  considering  each  of  the  tangent 
planes  to  the  torus,  the  critical  points  may  be  identified  as  those  at  which  tangent  planes  coincide 
with  height  planes  of  the  Morse  function.  For  example  the  torus,  as  we  expected  there  will  be 
critical  points  corresponding  to  the  maximum  and  minimum  (at  the  north  and  south  poles  of  the 
torus)  and  the  two  saddle  points  at  the  handle. 

There  is  further  geometric  interpretation  of  Morse  theory  for  a  2-manifold  by  correlating  the 
tangent  plane  of  the  surface  at  each  point  to  the  planes  defined  by  a  height  function.  Specifically, 
we  consider  classifying  critical  points  and  trivial  points  (non-critical  points).  Similar  to  the  method 
presented  above,  where  we  classify  points  based  on  their  shape,  we  consider  analyzing  the  local 
shape  of  the  surface,  by  looking  at  the  relationship  between  a  small  circular  neighborhood  of  each 
point  on  the  surface  and  the  height  planes  of  a  height  function. 

Morse  lemma  says  that  near  po  there  is  a  smooth  change  of  coordinates  under  which  the  resulting 
Taylor  series  of  the  Morse  function  h  near  Po  is  the  pure  quadratic  function. 

Theorem  7.7  Morse  functions  are  stable  and  dense  in  the  set  of  all  smooth  functions.  Equiva¬ 
lently,  any  smooth  function  can  be  converted  into  a  Morse  function  as  a  result  of  a  perturbation  as 
slight  as  desired. 

This  Morse’s  theorem  says  that  a  small,  smooth  perturbation  of  a  Morse  function  yields  another 
Morse  function.  The  density  means  that  there  is  a  Morse  function  arbitrarily  close  to  any  non-Morse 
function. 

7.2.5  Generalized  height  function 

The  height  function  in  the  direction  of  a  vector  v  (we  may  assume  v  €  S^,  i.e.  ||u||  =  1)  is  defined 
as  fiv  :  M  — >  R  such  that  hv{p)  =  p  •  v. 

The  level  sets  of  the  height  function  are  the  intersections  of  M  with  planes  orthogonal  to  v 
(considered  as  a  line)  as  pictured  in  Figure  7.4.  Denote  by  Pz  the  plane  at  height  z.  The  original 
object  (surface)  can  be  reconstructed  if  we  know  all  its  sections  by  these  parallel  planes  (i.e.  the 
surface  is  the  union  of  these  planes). 
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Figure  7.4:  Illustration  of  the  height  function. 

Clearly,  the  level  sets  of  the  height  function  may  have  isolated  points,  or  curves,  or  may  contain 
an  open  subset  of  the  plane.  Furthermore,  the  level  sets  may  be  connected  or  disconnected,  and 
the  curves  may  have  complicated  singularities  (e.g.  generalized  Monkey  saddle  surface). 

A  point  p  is  a  critical  point  of  h  if  and  only  if  g{p)  =  where  ±v  is  the  point  in 
corresponding  to  u  €  S^,  that  is,  the  inverse  image  g~^{±v)  consists  exactly  of  all  points  on  the 
surface  M  whose  tangent  planes  are  orthogonal  to  u.  Hence,  p  is  a  nondegenerate  critical  point  of 
h  if  and  only  if  it  is  a  regular  point  for  g.  Therefore,  rbv  is  a  regular  value  of  g  if  and  only  if  all  the 
critical  points  of  h  are  nondegenerate. 

Proposition  7.8  The  height  function  h  :M-^R  in  the  direction  o/  u  G  is  a  Morse  function  if 
and  only  if  the  corresponding  point  is  a  regular  value  for  the  Gauss  map  p  :  M  — >  RP^. 

Applying  Sard’s  theorem  to  the  map  p,  we  conclude  that  the  set  of  for  which  the  height  function 
h  is  not  a  Morse  function  has  measure  zero  in  RP^. 
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7.2.6  Height  function  and  immersion 

Defining  a  height  function  for  a  closed  2-manifold  and  examining  its  pre-image  for  various  intervals 
along  the  z-direction,  will  create  a  sequence  of  closed  contours  on  the  surface.  This  corresponds  to 
placing  the  closed  2-manifold  in  a  tank  and  slowly  immersing  it  in  liquid  up  to  various  heights  by 
adding  more  and  more  water  to  the  tank.  The  level  set  for  a  given  height  z  will  be  the  intersection 
of  the  surface  with  the  top  of  the  water  (see  Figure  7.5).  We  observe  that  as  the  surface  is  immersed 
in  the  water,  the  topology  of  the  level  sets  change,  i.e.,  the  number  of  components  of  the  level  set 
changes  for  various  heights.  For  example,  imagine  we  are  pouring  water  into  a  tank  with  the  surface 
shown  in  Figure  7.5.  As  we  first  pour  water  in  to  level  zq,  we  do  not  intersect  the  surface  and  the 
pre-image  of  our  height  function  will  be  empty  (it  will  have  no  contours).  When  the  water  first 
touches  the  surface  at  level  Zj,  the  topology  of  the  level  set  changes  and  the  pre-image  now  consists 
of  a  single  contour.  As  we  continue  pouring  water  into  the  tank,  the  topology  of  the  level  sets  will 
continue  to  change.  For  example,  when  the  water  level  first  reaches  a  “hole”  of  one  of  the  handles, 
the  topology  of  the  level  set  will  change  from  a  single  contoiu:  to  two.  Finally,  consider  when  we 
pour  in  the  last  of  the  water  and  the  level  set  changes  such  that  we  once  again  have  no  contours. 
This  analogy  of  immersing  a  surface  in  water  is  often  used  to  describe  the  process  of  finding  critical 
points  in  a  surface. 
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7.3.1  Reeb  graph 

An  interesting  concept  related  to  Morse  theory  and  very  useful  in  analyzing  a  surface  topology 
is  the  Reeb  graph.  The  latter  is  defined  as  a  quotient  space  M/~  with  an  equivalence  relation 
given  by  p  ~  q  if  and  only  if  h(p)  —  h{q)  and  p,  q  belonging  to  the  same  connected  component 
of  h~^(h{p)).  In  particular,  each  connected  component  is  represented  by  a  point  in  the  Reeb 
graph  as  illustrated  in  Figure  7.6.  The  left  figure  shows  a  torus  with  the  critical  points  of  its 
height  function  (Morse  function).  The  figure  in  the  middle  illustrates  the  geometric  features  of  the 
torus  represented  by  cross-sections  along  its  height.  The  right  figure  shows  the  topological  features 
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Figure  7.5:  Surface  immersed  in  water. 

captured  by  the  Reeb  graph.  By  taking  an  appropriate  number  of  cross-sections  and  a  smooth 
interpolation  in  between,  Shinagawa  et  aL  [67, 68]  proposed  a  Reeb  graph  based  approach  or  so- 
called  homotopy  model  for  object  reconstruction.  The  Reeb  graph  is  a  topological  representation  of 
an  object  (skeletal  structure),  and  has  storage  and  transmission  advantages  due  to  a  parsimonious 
data  representation. 

Mathematically,  a  quotient  space  M/~  =  {[p]  :  p  €  M}  is  a  set  of  equivalence  classes  with 
relation  and  where  [p]  =  {q  6  M  :  q  ^  p}  is  the  equivalence  class  of  p  E  M.  Intuitively,  M/^  is 
a  space  created  by  taking  the  space  M  and  gluing  p  to  any  q  that  satisfies  q  ~  p.  The  classes  [p]  are 
the  connected  components  for  the  Reeb  graph,  and  being  in  the  same  component  is  an  equivalence 
relation: 

q  p  h{q)  =  h{p)  and  p,  q  E  ConComp{/i“^(fe(p))}, 

where  ConComp{-}  denotes  the  connected  component.  In  a  Reeb  graph  representation,  each  con¬ 
nected  component  of  a  contour  (i.e.  h~^{z)  where  z  =  h{x,y^z))  corresponds  to  a  point. 


Figure  7.6:  Reeb  graph  representation  of  a  torus. 

7.4  Level  sets  around  Morse  points 

Let  /  :  M  ^  R  be  a  Morse  function  defined  on  a  compact  surface  M.  The  following  result  shows 
that  a  Morse  function  on  a  surface  may  determine  the  shape  of  the  surface. 

Proposition  7.9  If  f  :  M.  has  exactly  two  nondegenerate  singular  points,  then  M  is  diffeo- 
morphic  to  the  sphere  S^. 

Denote  by  Ma  the  sub-surface  of  M  consisting  of  all  points  at  which  h  takes  values  less  than  or 
equal  to  a  real  number  a 

Ma  =  {p  €  M  ;  /(p)  <  a} 

and  denote  by  La  the  set  of  points  where  the  value  of  h  is  exactly  o,  that  is  La  =  Note 

that  when  a  is  a  regular  value,  the  set  La  is  a  smooth  curve  of  M  and  it  is  the  boimdary  of  Mo- 

Proposition  7.10  Let  a  <  b  be  real  numbers  such  that  the  function  f  :  M.  —y  M.  has  no  critical 
value  in  the  interval  [a,  6] ,  then 
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Figure  7.7;  Illustration  of  Mq  and  La- 

(a)  The  level  curves  La  =  and  Lf,  =  are  diffeomorphic. 

(b)  The  subsurfaces  Ma  and  Mfc  are  diffeomorphic,  with  boundaries  La  and  Lb  respec¬ 
tively. 

Figure  7.8  shows  the  evolution  of  the  subsurface  Ma  as  the  parameter  a  changes.  If  a  < 
minpgM{/(p)}>  then  Mq  =  0.  And  as  we  increase  the  parameter  a,  the  subsurface  Ma  changes 
until  it  covers  the  entire  surface  M.  We  may  think  of  the  height  function  /  :  M  — ►  R  as  dipping  a 
doughnut  into  a  cup  of  chocolate  cream. 
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Figure  7.8:  Evolution  of  Mq  as  a  changes. 

7.4.1  Handle  decompositions 

Let  /  :  M  — »  E  be  a  Morse  function  defined  on  a  compact  surface  M.  Each  time  the  value  of  /  passes 
through  a  critical  value,  a  handle  appears  and  is  attached  to  the  previously  built-up  subsurface. 
The  index  of  the  handle  coincides  with  the  index  of  the  corresponding  critical  points,  that  is  the 
number  of  negative  eigenvalues  of  the  Hessian  matrix  of  /.  For  example,  let  Pi,P2fP3>P4  b®  fbe 
critical  points  of  a  height  function  /i :  T  — >  M  defined  on  a  torus  T,  and  denote  by  vi,V2,V3,  V4  their 
critical  values,  i.e.  h(pi)  =  Vi  is  the  ^-coordinate  of  each  point  p^.  We  further  assume  that  these 
critical  values  are  ordered  vi  <  V2  <  V3  <  V4,  that  is  ui  is  the  minimum  value,  V4  is  the  maximum 
value,  and  V3,  V4  are  the  saddle  values. 

To  track  the  topological  changes  of  the  surface  M,  we  look  at  how  Mq  changes  as  the  parameter 
a  increases.  In  the  case  of  a  torus,  we  start  from  a  value  less  than  vi,  that  is  for  a  <  we  have 
Mo  =  0.  As  soon  as  a  passes  Vi,  a  2-disk  (upright  bowl)  pops  out  and  we  have  Mo  =  D"^.  This 


7.5  Distance  function 


114 


2-disk  corresponds  to  the  minimum  critical  point  of  index  0  and  is  called  0~handle.  Similarly,  a 
1-handle  corresponds  to  a  saddle  point  and  a  2-handle  corresponds  to  a  maximum  critical  point. 

Consequently,  we  deduce  that  any  closed  surface  can  be  decomposed  into  a  union  of  a  finite 
number  of  0-,l-,  and  2-handles.  In  other  words,  any  closed  surface  admits  a  handle  decomposition. 

The  diagram  depicted  in  Figure  7.9  shows  the  sequence  of  steps  in  the  gradual  buildup  of  a 
torus,  starting  with  a  disk  (or  0-handle),  adding  two  consecutive  1-handles,  and  finally  completing 
the  torus  with  a  2-handle. 


Figure  7.9:  handle  decomposition. 


7.5  Distance  function 

The  concept  of  distance  is  central  to  topology,  with  the  actual  numeric  values  being  of  less  im¬ 
portance.  In  fact,  topologists  often  use  a  distance  function,  but  the  attributed  numerical  values 
have  only  secondary  meaning.  To  illustrate  this,  suppose  we  are  given  an  object  in  the  ordinary 
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three-dimensional  space,  and  a  point  outside  the  object,  and  the  question  is:  does  the  object  come 
arbitrarily  close  to  this  reference  point?.  This  may  be  stated  as:  is  the  point  a  boundary  point  of 
the  object?  “Arbitrarily  close”  means  that  if  one  imagines  a  ball  around  the  reference  point,  then 
the  ball  contains  some  points  belonging  to  the  object  no  matter  how  small  the  ball  is.  The  actual 
distances  between  the  points  belonging  to  the  object  and  the  reference  point  do  not  matter,  and 
there  just  have  to  be  arbitrarily  small  values  among  them. 

The  distance  function  is  a  function  which  has  nondegenerate  critical  points,  and  it  can  be  shown 
that  almost  all  distance  functions  are  Morse  functions.  In  fact,  for  fixed  v  €  M^,  we  niay  define  a 
distance  function  of  M  to  v  as  dt;  •  M  ►  R  such  that  dv{p)  =  ||p  —  v\\^. 

If  a  surface  M  is  given  in  parametric  form  r{x^y)  where  {x^y}  is  the  coordinate  system,  then 
the  distance  function  may  be  expressed  as  dx;(r(x,  j/))  =  ||r(rr,  y)  —  v\[^.  The  first  partial  derivatives 
are  given  by  dx  =  2rx  •  (r(x,  y)  —  q)  and  dy  =  2ry  •  (r(a;,  y)  —  q).  Hence  d  has  a  critical  point  at 
p  =  r{Xjy)  if  and  only  if  v  —  p  is  orthogonal  to  M  at  p,  i.e.  v  ~  r(x, y)  is  parallel  to  the  surface 
normal  N.  Thus  v  =  r(x,  y)  -h  a  iV. 

The  distance  function  from  the  origin  of  a  coordinate  system  is  given  by  d(x,  y)  =  ||r(x,  y)\\^  — 
x^  +  y^  +  u(x,  y)^.  Its  gradient  is  Vd(x,  y)  =  2[(x,  y)  +  u(x,  y) Vu(x,  y)],  and  its  Hessian  matrix  is 
V^u(x,y)  =  2[(1  +  ||Vn(x,y)p)72  +w(x,y)(V^w(x,y))],  where  h  is  the  2x2  identity  matrix.  The 
second  partial  derivatives  at  a  critical  point  can  be  easily  derived  as 

dxx  ~  2(t*/C  *  ^ X  ^ XX  *  (^(^5  y)  '*^))  ^  2(1*3.  *  Vx  UjTxx  * 

Hence,  the  Hessian  matrix  may  be  expressed  in  terms  of  the  first  and  second  fundamental  forms  as 
follows 

V^d  =  2{I-aJ[), 

where  I  and  I  are  the  first  and  second  fundamental  forms  respectively. 

A  degenerate  critical  point  of  the  distance  function  satisfies  det(V^d)  =  0  if  and  only  if 
det(V^w)  =  1/a^  =  «i«2j  where  ki  and  K2  are  the  principal  curvatures.  A  point  p  €  M  is 
therefore  a  degenerate  critical  point  of  the  distance  function  dy  if  and  only  if  v  is  a  focal  point  of 
(M,p).  In  addition,  the  Morse  index  of  a  nondegenerate  critical  point  of  the  distance  function  dy 
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is  equal  to  the  number  of  focal  points  of  (M,p)  which  lie  on  the  segment  from  p  to  v.  This  can  be 
shown  using  the  Hessian  matrix  since  the  number  of  its  negative  eigenvalues  is  equal  to  the 
number  of  eigenvalues  of  the  JT  (assuming  that  I  is  the  identity  matrix)  which  are  >  1/a. 

Without  loss  of  generality  we  choose  v  to  be  the  centroid  c  of  the  surface  M,  and  for  simplicity  we 
consider  the  centroid  to  be  the  origin  of  the  Euclidean  coordinate  system  as  pictured  in  Figure  7.10. 
Hence  the  distance  function  becomes 

d{p)  =  ||pf  =  +  + 

where  p  =  (re,  y,  2r).  Note  that  for  r  >  0,  the  level  sets  {p  G  M  :  d(p)  =  r^}  of  the  distance  function 
are  concentric  spheres  of  radii  r,  and  the  object  can  be  reconstructed  if  we  know  its  intersections 
with  these  concentric  spheres  (see  Figure  7.10). 


Figure  7.10:  Illustration  of  the  distance  function. 
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7.6  Connection  between  height  function  and  distance  function 

For  certain  purposes  such  as  terrain  image  reconstruction,  the  height  function  has  nicer  properties, 
while  for  others  the  distance  function  behaves  better  due  mainly  to  its  rotational  invariance.  There 
is,  however,  one  situation  when  these  two  functions  are  essentially  the  same.  Suppose  that  a  surface 
M  is  embedded  in  a  sphere  centered  at  the  origin  and  with  radius  R  (see  Figure  7.11),  then  the 
height  function  hy  and  the  distance  function  dy  differ  by  a  constant  and  therefore  have  the  same 
critical  points: 

dv{p)  =  llp-vf  =  +  IHI^  -2p-'u  =  {\\vf  +  R^)-2hv(p). 


Figure  7.11:  Embedding  of  a  3D  airplane  into  a  sphere. 

The  key  idea  behind  using  the  distance  function  is  to  track  the  changes  in  topology  as  we 
cross  a  surface  singularity.  In  the  first  step,  we  start  with  a  sphere  having  a  sufficiently  small 
radius,  and  centered  as  the  barycenter  of  the  underlying  surface,  then  we  evolve  this  sphere  by 
increasing  its  radius  so  that  we  will  have  a  set  of  concentric  spheres  covering  the  entire  surface.  As 
we  cross  a  surface  singularity,  a  topological  change  of  the  level  curves  will  take  place  as  illustrated 
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Figure  7.12:  Distance  function  defined  on  a  torus. 


Figure  7.13:  Distance  function  defined  on  a  torus  (cont.). 

in  Figures  7.12,  7.13,  and  7.14.  The  level  curves  for  two  3D  real  data  objects  are  pictured  in 
Figure  7.15.  In  other  words,  some  level  curves  may  split  or  merge.  We  are  essentially  interested 
in  these  changes  for  topological  modelling  purposes  which  in  turn  may  be  explained  by  applying 
Morse  theory  to  the  distance  function. 


Chapter  8 


Conclusions  and  Future  Research 

This  thesis  has  presented  computational  algorithms  for  variational  image  denoising,  topological 
modeling,  three-dimensional  object  recognition,  and  geometric  matching.  We  have  demonstrated 
the  use  of  these  algorithms  through  a  variety  of  imaging  and  computer  vision  applications  including 
image  filtering,  singularity  extraction  and  evolution,  topological  modeling  of  illuminated  surfaces, 
geodesic  matching  of  triangle  meshes,  and  distance  function-based  object  recognition.  The  geomet¬ 
ric/topological  algorithms  are  tailored  for  the  discrete  representation  of  surfaces  as  triangle  meshes. 
We  have  demonstrated  the  effectiveness  of  the  proposed  methods  through  numerical  simulations 
with  synthetic  and  real  data  in  2D  and  3D  computer  imagery. 

In  the  next  Section,  the  contributions  made  in  each  of  the  previous  chapters  and  the  concluding 
results  drawn  firom  the  associated  research  work  are  presented.  Suggestions  for  future  research 
directions  related  to  this  thesis  are  provided  in  Section  8.2. 

8.1  Contributions  of  the  thesis 

8.1.1  Robust  and  efficient  variational  filters  for  image  denoising 

Using  the  theoretical  fundamentals  of  robust  statistics,  a  variational  filter  referred  to  as  a  Huber 
gradient  descent  flow  was  proposed  in  Chapter  3.  It  is  a  result  of  optimizing  a  Huber  functional 
subject  to  some  noise  constraints,  and  it  takes  a  hybrid  form  of  a  total  variation  diffusion  for  large 
gradient  magnitudes  and  of  a  linear  diffusion  for  small  gradient  magnitudes.  Using  the  gained 
insight,  and  as  a  further  extension,  we  proposed  an  information-theoretic  gradient  descent  flow 
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which  is  a  result  of  minimizing  a  functional  that  is  a  hybrid  between  a  negentropy  variational 
integral  and  a  total  variation.  Illustrating  experimental  results  demonstrated  a  much  improved 
performance  of  the  approach  in  the  presence  of  Gaussian  and  heavy-tailed  (impulsive)  noise. 

8.1.2  A  topological  variational  model  for  image  singularities 

Image  singularities  are  prominent  landmarks  and  their  detection,  recognition,  and  classification 
is  a  crucial  step  in  image  processing  and  computer  vision.  Such  singularities  carry  important 
information  for  further  operations,  such  as  image  registration,  shape  analysis,  motion  estimation, 
and  object  recognition.  In  Chapter  4,  we  proposed  a  topological  gradient  descent  flow  for  image 
singularities.  The  approach  is  expressed  in  the  higher  order  variational  framework  as  a  minimizer  of 
a  variational  integral  involving  the  gradient  and  the  Hessian  matrix  of  the  height  function  defined  on 
a  manifold.  We  demonstrated  through  numerical  simulations  the  power  of  the  proposed  technique 
in  preserving  image  singularities. 

8.1.3  Topological  modeling  of  illuminated  surfaces  using  Reeb  graph 

In  Chapter  5  we  introduced  a  reliable  and  efficient  feature  based  object  representation  for  topological 
modeling  of  three-dimensional  illuminated  surfaces.  The  proposed  approach  encodes  an  object 
into  the  Reeb  graph  concept  from  computational  topology.  This  skeletal  structure  is  based  on  a 
generalized  height  function  in  the  light  direction  defined  on  an  illuminated  surface.  The  topological 
properties  of  the  proposed  representation  were  analyzed  in  the  Morse  theoretic  framework,  and  its 
close  relationship  to  the  shading  problem  was  also  highlighted.  Some  numerical  simulations  with 
S3Tithetic  and  real  3D  data  were  provided  to  demonstrate  the  potential  of  object  singularities  in 
topological  modeling. 

8.1.4  Geodesic  object  representation  and  recognition 

In  Chapter  6  we  proposed  a  shape  signature  that  captures  the  intrinsic  geometric  structure  of 
3D  objects.  The  primary  motivation  of  the  proposed  approach  is  to  encode  a  3D  shape  into  a 
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one-dimensional  geodesic  distribution  function.  This  compact  and  computationally  simple  repre¬ 
sentation  is  based  on  a  global  geodesic  distance  defined  on  an  object  surface,  and  takes  the  form  of 
a  kernel  density  estimate.  To  gain  further  insight  into  the  geodesic  shape  distribution  and  its  prac¬ 
ticality  in  3D  computer  imagery,  some  numerical  experiments  were  provided  to  demonstrate  the 
potential  and  the  much  improved  performance  of  the  proposed  methodology  in  3D  object  matching. 
This  was  carried  out  using  an  information-theoretic  measure  of  dissimilarity  between  probabilistic 
shape  distributions. 

8.1.5  Distaoice  function-based  object  recognition 

In  Chapter  7  we  introduced  a  topological  approach  to  object  recognition  using  a  distance  function. 
Similarly  to  the  height  function  strategy  which  consists  of  reconstructing  surface  from  its  cross- 
sections,  the  key  idea  behind  using  a  distance  function  is  that  a  surface  may  be  reconstructed  from 
its  intersections  with  concentric  spheres  centered  at  the  centroid  of  the  underlying  surface.  The 
topological  changes  in  the  surface  occur  as  we  increase  the  value  of  the  sphere  radius.  At  singular 
points,  the  level  curves  of  the  distance  function  may  split  or  merge  which  indicate  topological 
changes.  We  also  show  that  when  a  surface  is  embedded  in  a  sphere,  the  height  function  and  the 
distance  function  are  equivalent  in  a  Morse-theoretic  setting,  that  is  both  functions  have  the  same 
singularities. 

8.2  Future  resecurch  directions 

Several  interesting  research  directions  motivated  by  this  thesis  are  discussed  next.  In  addition  to 
designing  new  methodologies  for  image  denoising  and  segmentation,  we  intend  to  accomplish  the 
following  projects  in  the  near  future: 

8.2.1  Attributed  Reeb  graph  matching,  indexing,  and  retrieval 

Recently  we  have  been  working  on  the  representation,  matching,  indexing  and  retrieval  in  a  database 
of  3D  objects  based  on  the  topological  and  geometric  information.  Building  a  database  requires 
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collecting  3D  models,  computing  their  Reeb  graph  representations,  and  indexing  in  the  data  base 
based  on  an  abstracted  information  given  by  their  Reeb  graphs.  An  appropriate  and  efficient  rep¬ 
resentation  of  the  Reeb  graph  is  the  attributed  Reeb  graph  that  represents  topology  and  geometry 
in  a  compact  representation,  where  vertices  and  edges  have  geometric  attributes.  In  other  words, 
we  associate  to  the  graph  as  much  geometric  information  as  possible  that  will  be  attached  to  the 
graph  for  further  tasks  such  as  matching,  indexing  and  retrieval. 

8.2.2  Entropic  minimum  spanning  Reeb  trees  for  terrain  image  analysis 

The  vertices  of  a  Reeb  tree  can  be  characterized  using  the  minimum  spanning  tree  (MST)  which 
aims  to  quantify  spacial  dot  patterns  by  revealing  hidden  nearest-neighbor  correlations:  The  MST 
representation  is  naturally  translation  and  rotation  invariant,  and  therefore  constitutes  a  good  can¬ 
didate  for  geo-registration  and  other  image  registration  applications.  The  Jensen-Renyi  divergence 
may  be  used  as  a  robust  dissimilarity  measure  between  the  Morse  features  of  the  target  and  the 
reference  images. 

8.2.3  Divergence  measures  and  information  geometry 

There  are  many  possibilities  for  extending  the  Jensen-Renyi  divergence  using  more  generalized 
entropy  measures,  and  borrowing  concepts  from  information  geometry  in  order  to  fully  employ  the 
mathematical  machinery  of  differential  geometry  and  topology.  Information  geometry  is  the  branch 
of  probability  theory  dedicated  to  provide  families  of  probability  distributions  with  differential 
geometrical  structures.  One  then  uses  the  tools  of  differential  geometry  in  order  to  have  a  clear 
and  intuitive  picture  of  a  family  of  probability  distributions  which  form  a  differentiable  manifold. 
Information  geometry  elucidates  the  geometric  structure  of  such  a  manifold. 

Theoretically  we  hope  to  develop  more  rigorous  analysis  for  the  Jensen-Renyi  divergence.  It 
is  also  worthwhile  to  combine  other  computational  techniques  with  our  approach.  Furthermore, 
we  are  planning  to  apply  this  divergence  measure  to  other  imaging  applications  including  DNA 
segmentation,  microarray  images  and  independent  component  analysis. 
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The  first  variation  of  the  functional  J^{u)  =  Jq  F(|  Vu|)  dx  in  the  direction  of  v  is  given  by 

. f  dx. 


6J^{u]  v)  =  +  eu) 

de 


■n-, 


The  following  identity 


yields 


F'(|Vn|) 


Vu  •  Vv  dx 


Jn  |V«|  |V«| 

Using  the  divergence  theorem  for  a  vector  field  w 


div(nVu)  =  div(V«)w  +  Vu  •  Vv, 


F'(|V«|), 


Vu  I  dx 


I  div(«;)  dx  =  I  w  -uds 

Ju  Jdn 


where  i/  is  the  outward  unit  normal  vector  (field)  on  dUl  (the  boundary  of  fi)  and  ds  is  an  area 


element.  Therefore 


Jst  V  |V»I  J  Jat 


h  \  |V«I  J  Jm  |V« 

If  we  assume  homogenuous  Neumann  boundary  conditions 

—  du 

vu  •  1/  =  wj/  =  —  =  0, 
ou 

then  the  first  variation  of  J'  is  reduced  to 


5T{u-v)  =  - 


which  concludes  the  proof. 
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A  numerical  implementation  of  the  partial  differential  equation  given  by  Eq.  13  is  performed  using 
an  explicit  scheme  in  time  and  location  as  follows.  Let  be  the  approximation  of  u{x,  y,  t)  on  a 
grid  {iAx,jAy,nAt).  For  simplicity  we  assume  that  Ax  =  Ay  =  h.  Denote  by 


q,n  _  „n 

_ I 


Dlu  =  ±- 


the  matrices  of  column  differences  and  row  differences  respectively  (i.e.  backward  and  forward 
differences). 

Similarly,  the  central  differences  are  given  by 


The  operator  div(g'(|  V«|)  V«)  on  the  right  hand  side  of  Eq.  13  is  discretized  using  an  upwind  scheme 
as  follows 


div^5(|Vti|)V«^  = 


where  minmod  is  a  function  that  returns  the  argument  with  smallest  absolute  value  when  all  the 
arguments  are  of  the  same  sign  and  zero  otherwise.  The  minmod  function  is  a  limiter  whose  goal 
is  to  prevent  oscillations  while  maintaining  the  order  of  accuracy  of  the  method,  and  it  is  defined 
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as 

j/’  1.N  /sign(a) +  sign(6)\  . 

minmod(a,  6)  =  ( - - - 1  min(|a|,  |6|) 

=  min(max(o,  0),  max(6, 0))  +  max(min(a,  0),  min(6, 0)) 

0  if  o6  <  0 

a  if  |a|  <  |fe|  and  ab  >  0 
b  if  |al  >  |6|  and  ab  >  0. 

Figure  (a)  below  depicts  the  gradient  vector  field  of  the  minmod  function  (the  background  color 
indicates  the  value  of  the  minmod  function),  and  Figure  (b)  illustrates  its  contours. 


(a)  3D  plot  of  the  minmod  function 


(b)  level  curves 


'appendix  C_ 

Appendix  C 


Surface  Curvatures 

Let  p  6  M.  Let  v  €  TpM  =  span{rx,rj,}  (we  assume  v  a  unit  vector).  By  definition  v  €  TpM  if 
there  exists  a  curve  c  :  (— c,  e)  — »  M  for  some  e  >  0  such  that  c(0)  =  p  and  c'{0)  =  v.  In  other  words, 
the  tangent  space  is  the  set  of  vectors  orthogonal  to  the  surface  normal.  Thus,  the  tangent  plane 
at  Po  is  the  set  of  points  p  such  that  JV(p)  •  (p  —  Po).  Hence  for  the  Monge  patch,  the  equation  of 
the  tangent  plane  at  a  point  (xp,  yo,  u{xo,  yo)  is  given  by  2:  =  u{xo,  yo)  +  Vu{xo,  y(i)-{x  —  XQ,y-  yo). 

Let  /  :  M  — >  M  be  a  smooth  function.  The  directional  derivative  of  /  at  p  in  the  direction  of  v 
is  given  by  Vvf  =  {f  °  c)'(O),  where  c  :  (— e,  e)  — >  M  for  some  e  >  0  is  a  curve  such  that  c(0)  =  p 
and  (0)  =  V. 

If  /  :  M  — >  R^,  then  the  differential  of  /  at  p  is  the  map  df  :  TpM  — >  R^  such  that  df{p)  =  V^/ 
for  all  v  €  TpM. 

If  /  :  M  — M  is  a  smooth  function  between  two  manifolds,  then  the  differential  of  /  at  p  is  the 
linear  map  df  :  TpM  — >  Ty^pjM.  In  particular,  let  v  E  TpM.  Then  G  TpM.  Indeed,  since 

<N,N  >=  1,  it  follows  that  0  =  Vv  <  iV,  iV  >=  2  <  V^iV,  N  >.  Thus,  VvN  is  orthogonal  to 

N. 

The  first  fundamental  form  I :  TpM  x TpM  — ♦  R  is  a  bilinear  form  such  that  I(v,  w)  —<  v,  w  >, 
and  its  matrix  with  respect  to  the  orthogonal  basis  {rx,rj,}  6  TpM  is  given  by 

I{rx,rx)  I{rx,ry) 

I{ry,rx) 
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The  matrix  I  is  also  called  the  metric  or  metric  tensor,  and  it  is  the  analogous  of  the  speed  in 
the  case  of  space  curves.  This  matrix  determines  all  the  intrinsic  properties  of  the  surface.  These 
properties  depend  on  the  surface  and  do  not  depend  on  its  embedding  in  space.  Furthermore,  the 
matrix  I  is  invariant  to  rotation  of  the  surface  in  space  because  it  is  defined  in  terms  of  inner 
products  that  are  rotation  invariant. 

The  Weingarten  map  (also  called  the  shape  operator)  is  a  linear  map  W  :  TpM.  — >  TpM 
such  that  W{v)  =  —  V^iV.  In  terms  of  the  basis  {rx,ry)  €  TpM,  we  have  W{rx)  =  —Nx  and 
W{ry)  =  —Ny.  These  equations  are  called  Weingarten  equations  and  they  express  the  derivatives 
of  the  normal  to  a  surface  using  derivatives  of  the  position  vector. 

The  second  fundamental  form  E  :  IpM  x  TpM  ^  R  is  a  bilinear  form  such  that  E(v,w)  = 
W(v)  ■  w  =  —  <  VvNyW  >,  and  its  matrix  with  respect  to  the  orthogonal  basis  {ra;,ry}  €  TpM 
is  given  by 

E{rx,rx)  E{rx,ry)  \  _  I  -Nx-rx  -Nx-Vy 

E{ry,rx)  E(ry,ry)  j  y  -Ny-Vx  ~Ny’Vy 

The  matrix  M  is  also  invariant  under  rotation  of  the  surface. 

The  third  fundamental  form  M  :  TpM  x  TpM  — >  R  is  a  bilinear  form  such  that  E[(v,  w)  =  W {v)  • 
W(tn)  =<  Vt;A^,  VtylV  >,  and  its  matrix  with  respect  to  the  orthogonal  basis  {ra;,ry}  €  TpM  is 
given  by 

M{rx,rx)  M{rx,ry) 

Er{ry,rx)  E[{ry,ry) 

The  third  fundeimental  form  is  given  in  terms  of  the  first  and  second  forms  hy  M  —  2HE  +  KI  =  0, 
where  H  and  K  are  the  mean  and  Gaussian  curvatures  respectively. 

The  Gaussian  curvature  is  given  by  the  determinant  of  W.  The  normal  curvature  of  M  in  the 
direction  v  6  TpM  is  k{v)  =  W{v)-v.  Since  v  G  TpM  =  span{ra;,  r^},  it  follows  that  v  =  arx+bry. 
Therefore  the  normal  curvature  in  the  direction  v  can  be  expressed  as 


k{v)  =  W{v)-v  = 


ia?  +  2mab  +  n6^ 
ea^  +  2fah  +  56^  ’ 


where  e,  /  and  g  axe  the  coefficients  of  the  first  fundamental  form,  and  £,  m  and  n  are  the  coefficients 
of  the  second  fundamental  form.  The  maximum  and  minimum  values  of  the  normal  curvature  at  a 
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point  on  a  regular  surface  are  called  the  principal  curvatures  ki  and 

Rodrigues  curvature  formula  is  given  by  dN  +  Kidr  =  0,  where  Ki  are  the  principal  curvatures. 
The  harmonic  curvature  of  the  principal  curvatures  is  defined  as 


The  extended  Gaussian  image  is  the  reciprocal  of  the  Gaussian  curvature,  that  is 


C{p)  - 


i 


Implicit  surface  curvatures 

Let  U  :  r2  C  — ♦  K  be  a  smooth  function.  *A.ii  implicit  surface  is  defined  as  a  level  set  of  the 

function  U,  i.e.  U{x,y.z)  —  0  (for  instance).  The  principal  curvatures  and  the  principal  directions 
of  the  level  surface  satisfy  the  following  ecpiation 

where  /a  is  a  3  x  3  identity  matrix. 

Offset  surface 

Let  M  C  be  a  regular  surface.  The  surface  parallel  to  M  at  a  distance  t  >  0  is  the  set 
M  =  {q d{q,M)  ^  t) 

The  offset  surface  patch  r  at  distance  t  to  a  surfacxi  M  parametrized  by  r  is  defined  as 

r{x,  y)  —  r{x.  y)  +  t  N[x,  y),  (1) 

where  t  is  the  distance  of  the  parallel  surface'  from  the  original  one,  and  N  is  the  unit  surface 
normal  of  r. 


